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Abstract 

The  terms  renormalization  and  renormalization  group  are  explained  by  reference 
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is  then  discussed;  firat  as  a  comprehensive  review  and  second  concentrating  on  the 
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Part  I.  Introduction 


Part  I 

Introduction 

1  Preliminary 

The  term  “RG”  (or  “RNG”)  has  now  become  familiar  to  the  turbulence  modeling 
and  computational  fluid  dynamics  (CFD)  community.  However,  understanding  just 
what  RG  is  actually  about  is  often  viewed  as  an  obscure  mathematical  problem. 
Accordingly,  we  feel  that  there  is  a  need  for  an  appropriate  article  for  this  diverse 
community.  The  present  review  article,  which  surveys  many  aspects  of  applications 
of  RG  to  turbulence,  is  intended  to  fill  this  gap. 

Turbulence  in  fluids  has  presented  a  notoriously  difficult  problem  for  more  than  a 
century.  As  discussed  in  standard  text  books  on  fluid  motion,  the  major  difficulty  is 
the  simultaneous  existence  of  many  different  space  and  time  scales  at  high  Reynolds 
number.  All  these  scales  are  of  equal  importance,  hence  the  theoretical  difficulty. 

Now,  RG  was  a  tool  originally  developed  for  dealing  with  such  problems.  Indeed, 
the  work  of  Wilson  in  developing  the  RG  method  for  critical  phenomena  led  to  his 
being  awarded  the  Nobel  prize  in  Physics  in  1982.  In  The  Nobel  Prize  Winners: 
Physics  (Magill  1989),  we  read: 

“Wilson  clarified  a  natural  phenomenon  that  has  puzzled  scientists  throughout  the 

centuries:  the  behavior  of  substances  at  critical  points,  or  phase  transitions.  With  his 
renormalization  group  theory,  Wilson  divided  this  seemingly  insoluble  problem  into 
a  number  of  small  ones.” 

1.1  General  remarks 


As  we  have  said,  at  high  Reynolds  numbers  (Re),  a  wide  range  of  length  and  time 
scales  in  turbulent  flows  become  excited.  Direct  numerical  simulations  (DNS)  -  which 
are  three-dimensional,  time-dependent  numerical  solutions  in  which  all  significant 
scales  of  motion  are  computed  without  any  modeling  -  are  currently  restricted  to  low 
turbulence  Reynolds  numbers  and  simple  geometries  (Reynolds  1990).  Some  recent 
DNS  of  homogeneous  and  isotropic  turbulence  at  5123  resolutions  are  Chen  et  al. 
(1993),  Wang  et  al.  (1996),  and  Yeung  and  Zhou  (1997). 

Although  DNS  can  be  extremely  useful  in  many  areas  related  to  the  study  of 
turbulence  physics  and  the  assessment  of  theories  (ibid),  virtually  all  scientific  and 
engineering  calculations  of  nontrivial  turbulent  flows,  at  high  Re,  are  based  on  some 
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type  of  modeling  (Speziale  1991).  Most  modeling  efforts  have  been  focused  on  subgrid 

scale  models  for  large-eddy  simulations  and  Reynolds  stress  models  -  the  two  forms 
of  modeling  that  will  be  discussed  in  this  article. 

Large-eddy  simulations  (LES)  are  a  logical  ‘modeling  extension’  of  DNS.  They 
are  based  on  the  observation  that  the  small  scales  are  more  universal  in  character 
than  the  large,  energy  containing  scales  of  motion.  In  LES,  the  three-dimensional 
time-dependent  motion  of  these  large  scales  are  computed  directly  while  the  effects  of 
the  small  scales  on  the  large  scales  are  modeled.  This  procedure  leads  to  the  so  called 
‘subgrid  scale  modeling’  problem;  see  Rogallo  and  Moin  (1984),  (Reynolds  1990), 
Lesieur  and  Metail  (1996),  and  Moin  (1997).  The  traditional  concept  of  eddy  viscosity 
(Smagorinsky,  1963)  is  that  the  subgrid  turbulence  simply  enhances  the  viscous  term 
in  the  momentum  equation,  and  that  the  value  of  the  enhanced  eddy  viscosity  can 
be  determined  by  equating  the  resulting  increase  in  the  energy  dissipation  rate  to  the 
rate  at  which  energy  is  transferred  from  resolved  to  subgrid  scales.  LES  is  useful  in 
the  study  of  turbulence  physics  at  high  Reynolds  numbers  that  are  unattainable  by 
DNS. 

LES  is  also  intended  to  be  useful  in  the  development  of  turbulence  models  for 
the  prediction  of  the  complex  flows  of  technical  interest  where  simpler  modeling  ap¬ 
proaches  fail  (Reynolds  1990).  However,  all  applications  of  LES  have  been  in  simple 
geometries  where  Van  Driest  damping  could  be  used  -  an  empirical  approach  that 
generally  does  not  work  well  when  there  is  flow  separation  (Speziale  1991). 

We  note  that  LES  and  DNS  are  recent  developments,  being  direct  products  of 
the  rapid  increase  in  computer  capabilities.  The  Reynolds  stress  modeling  concept, 
however,  was  introduced  by  Osborne  Reynolds  as  long  ago  as  1895.  Reynolds  stress 
models  remain  the  most  important  tool  for  the  highly  repetitive  engineering  calcula¬ 
tions  associated  with  design  (Reynolds  1990).  The  goal  of  Reynolds  stress  modeling 
is  to  provide  reliable  information  about  first  and  second  one-point  moments  (e.g.,  the 
mean  velocity,  mean  pressure,  and  turbulence  intensity),  which  is  usually  all  that  is 
needed  for  design  purposes  (Speziale  1991).  In  principle,  this  approach  attempts  to 
model  all  of  the  fluctuating  scales. 

Recently,  the  methodology  of  RG  has  attracted  considerable  attention  as  a  sys¬ 
tematic  approach  to  subgrid  scale  modeling  in  LES  and  Reynolds  stress  modeling.  In 
principle,  the  RG  removal  of  only  the  smallest  scales  generates  a  subgrid  scale  model 
while  Reynolds  stress  models  are  obtained  in  the  limit  when  all  fluctuating  scales  are 
removed.  This  point  of  view  is  due  to  Yakhot  and  Orszag  (1986). 

1.2  Overview 

Because  we  are  addressing  such  a  diverse  community,  some  innovation  in  the  structure 
of  the  review  is  indicated.  Part  II  is  intended  for  those  readers  who  wish  to  obtain 
a  basic  knowledge  of  “RG”  and  what  are  the  main  results  from  this  approach  to 
studying  turbulence;  and  in  that  sense  it  is  complete  within  itself. 
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Part  III  is  for  those  readers  who  would  like  to  go  more  deeply  into  RG  theory. 
Here  we  present  the  detailed  mathematical  procedures  of  two  RG  approaches.  In 
particular,  we  clarify  their  major  approximations  for  the  benefit  of  the  reader  who  is 
willing  to  pursue  this  line  of  research.  In  this  part,  we  have  included  several  related 
pieces  of  work  and  have  also  attempted  to  compile  a  comprehensive  reference  list. 

Some  turbulence  modelers  or  CFD  workers  may  wish  to  skip  Part  III  and  proceed 
directly  to  Part  IV.  Here  we  present  recent  efforts  to  develop  turbulent  Reynolds 
stress  models  based  on  RG. 
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Part  II 

General  description  of  RG 

2  What  Is  Renormalization  Group? 

In  this  section  we  shall  explain  something  of  the  origins  and  meaning  of  RG.  In  the 
process  we  shall  define  terms  and  introduce  concepts  which  will  be  needed  when  we 
consider  the  question  of  how  RG  can  be  applied  to  fluid  turbulence. 

We  should  begin  by  noting  the  difference  between  the  concept  of  renormalization 
and  renormalization  group.  The  renormalization  procedure  was  introduced  originally 
as  a  method  of  removing  divergences  in  quantum  field  theory.  The  renormaliza¬ 
tion  group  theory,  on  the  other  hand,  was  used  for  improving  perturbation  theory 
(Gell-Mann  and  Low,  1953)  by  exploiting  the  non-uniqueness  in  the  renormalization 
procedure  (Stuechelberg  and  Petermann,  1953).  We  can  illustrate  this  idea  by  con¬ 
sidering  a  problem  which  may  be  discussed  entirely  in  terms  of  classical  physics.  This 
is  the  electron  gas,  which  consists  of  a  gas  of  electrons  in  a  uniform  background  of 
positive  charge,  so  that  the  system  is,  overall,  electrically  neutral.  Such  a  gas  can  be 
a  model  for  the  conduction  electrons  in  a  metal,  for  a  plasma,  or  for  an  electrolyte. 
This  allows  us  to  introduce  the  general  idea  of  renormalization,  and  is  the  subject  of 
Section  2.1.  In  order  to  introduce  the  more  specific  topic  of  renormalization  group, 
we  consider  as  an  example  the  phase  transition  from  para-  to  ferromagnetism.  This  is 
dealt  with  in  Section  2.2.  We  will  review  the  applications  of  RG  to  critical  phenomena 
in  Sections  2.3. 

2.1  Renormalization 

It  is  a  cardinal  assumption  of  physics  that  the  macroscopic  properties  of  a  material, 
such  as  its  dielectric  constant  or  elastic  modulus,  will  not  in  general  depend  on  the 
size  of  the  sample  of  material  under  consideration.  So,  for  a  microscopic  theory,  this 
implies  that  we  must  be  able  to  take  limits,  corresponding  to  the  distance  between 
constituent  microscopic  particles  shrinking  to  zero,  or  to  the  overall  size  of  the  system 
going  to  infinity,  without  upsetting  our  theory.  However,  if  we  take  our  example  of 
an  electron  gas,  it  is  easily  seen  that  problems  may  arise  in  both  these  limits. 

Let  us  suppose,  for  instance,  that  we  wished  to  work  out  the  total  electrostatic 
potential  energy  of  all  the  electrons  in  the  system  and  that  we  did  that  naively  by 
adding  up  the  contributions  from  individual  electrons,  each  electron  contributing  the 
usual  Coulomb  T/r’.  It  is  obvious  that  a  continuum  limit  would  require  r  —>  0,  with 
each  potential  diverging  in  that  limit.  Less  obviously,  the  limit  of  infinite  system  size 
also  leads  to  a  divergence.  Suppose  we  take  the  test  charge  to  be  at  the  origin  and 
assume  that  the  surrounding  charges  have  a  uniform  density  p,  say.  Then,  as  the 
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number  of  electrons  in  a  spherical  shell  between  r  and  r  +  dr  goes  as  pr2dr,  it  is  easily 
seen  that  this  overpowers  the  Coulomb  potential  to  give  a  total  potential  which  varies 
as  2 irpeR2,  where  e  is  the  electronic  charge  and  R  is  the  system  size. 

The  first  of  these  divergences  is  usually  referred  to  as  being  ‘ultra-violet’ ,  because 
short  distances  correspond  to  high  spatial  frequencies.  Similarly,  the  long-range  di¬ 
vergence  in  the  second  case,  is  referred  to  as  being  an  ‘infra-red  divergence’.  From  our 
point  of  view,  the  second  case  is  the  more  interesting.  In  a  charge  neutral  plasma,  the 
energy  at  a  point  does  not  diverge  with  system  size  because  a  cloud  of  ions  around 
each  individual  electron  acts  to  screen  the  long-range  Coulomb  potential  and  turns 
it  into  an  exponentially  decaying  potential  (for  stationary  particles).  Theoretically, 
this  effect  Bhows  up  when  we  take  the  interactions  between  individual  electrons  into 
account.  In  perturbation  theory  it  would  correspond  to  summing  a  certain  class  of 
terms  to  all  orders.  In  fact,  the  potential  at  lowest  nontrivial  order,  the  potential 
becomes  eexp[— r//D]/r,  where  lD  is  the  Debye  length  and  depends  on  the  density  of 
electrons. 

This  is  known  as  a  screened  potential.  But,  an  alternative  interpretation  is  to  say 
that  we  replace  the  electron  charge  e  by  the  renormalized  electron  charge  eexp[— r//D], 
This  replacement  of  a  ‘bare’  particle  by  a  ‘dressed’  particle,  in  order  to  take  account 
of  the  effect  of  interactions,  is  an  important  technique  in  many-body  physics  and 
nowadays  is  widely  referred  to  as  renormalization.  It  should  be  noted  that  a  charac¬ 
teristic  feature  of  the  process  is  that  the  renormalized  properties  (as  in  our  example) 
are  scale-dependent. 

This  discussion  serves  to  introduce  the  idea  of  renormalization.  Our  next  step  is 
to  consider  what  is  meant  by  renormalization  group ,  and  it  turns  out  that  this  is  just 
one  way  of  achieving  renormalization. 


2.2  Renormalization  group 

RG  was  based  on  the  obvious  equivalence  between  taking  the  continuum  limit  and 
making  a  scaling  transformation  of  the  basic  variables  such  as  position  or  momentum. 
However,  this  theory  was  regarded  as  rather  opaque  in  terms  of  its  underlying  physics, 
and  from  our  point  of  view,  it  is  better  to  consider  its  more  modern  incarnation  in  the 
theory  of  critical  phenomena,  in  which  it  is  seen  as  a  way  of  dealing  with  problems 
involving  a  large  range  of  length  or  time  scales.  Its  possible  relevance  to  turbulence 
might  now  begin  to  emerge;  but  we  shall  first  consider  another  simple  physical  system: 
the  lattice-spin  model  of  a  ferromagnet. 

Magnetism  arises  because  spins  at  different  lattice  sites  tend  to  become  aligned. 
This  tendency  is  opposed  by  thermal  effects,  which  promote  random  orientation  of 
the  spins.  At  any  fixed  temperature,  alignments  occur  on  length  scales  ranging  from 

the  lattice  spacing  L0  up  to  some  correlation  length  £.  Take  a  temperature  T0  which 
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is  sufficiently  high  that  £(T0)  is  of  order  L0.  It  is  evident  that  £  will  increase  as  T 
decreases  below  T0.  In  fact,  £  actually  becomes  infinite  as  T  approaches  the  Curie 
temperature,  which  defines  the  critical  point  for  this  problem.  At  this  critical  point, 
fluctuations  on  all  scales  from  L0  up  to  the  size  of  the  material  specimen  can  occur, 
and  therefore  a  net  overall  magnetization  can  appear. 

The  theoretical  problem  is  to  calculate  the  partition  function  and  thence  the 
thermodynamic  properties  of  the  material.  The  partition  function  is  a  sum  over  all 
possible  spin  configurations;  although  this  sum  is  easily  evaluated  if  the  spins  at 
different  sites  are  independent,  coupling  between  spins  makes  this  problem  difficult. 
The  fact  that  all  length-scales  are  (in  principle)  equally  important  makes  it  difficult 
to  know  what  to  do. 

The  problem  is  eased  to  some  extent  by  considering  models,  such  as  the  Ising 
model,  in  which  spin  vectors  are  Boolean  in  character  and  are  either  ‘up’  or  ‘down’, 
with  no  intermediate  states  permitted.  Another  helpful  feature  of  this  Ising  model 
is  that  the  spins  only  have  nearest-neighbour  interactions.  Then  the  application  of 
RG  to  this  problem  can  be  understood  as  a  form  of  coarse-graining,  in  the  following 
sense.  The  version  of  the  RG  method,  called  the  block-spin  transformation,  consists 
in  breaking  down  an  intractable  problem  with  multiple  scales  of  length  into  a  sequence 
of  smaller  problems,  each  of  which  is  confined  to  a  single  length  scale.  This  procedure 
consists  of  three  steps.  First,  the  lattice  is  divided  into  blocks  of  a  few  spins  each 
(A2,  say,  for  a  two-dimensional  lattice).  Then  each  block  is  replaced  by  a  single  spin 
whose  value  is  the  average  of  all  the  spins  in  the  block;  there  the  average  is  determined 
by  the  majority  rule.  Therefore,  a  new  lattice  is  created  with  N  times  the  original 
spacing  and  1  /N  the  density  of  the  spins.  Finally,  the  original  scale  is  restored  by 
reducing  all  dimensions  by  a  factor  of  N. 

We  now  provide  a  brief  description  on  the  RG  transformation  with  the  block 
spin  method.  Our  starting  point  is  the  Hamiltonian  H0,  associated  with  two  spins 
separated  by  a  distance  L0  {i.e.  the  lattice  spacing).  Then  we  calculate  an  effective 
Hamiltonian  Hi,  associated  with  regions  of  size  2Lo,  by  averaging  out  the  effects  of 
scales  Lq.  Next  we  calculate  Hi,  associated  with  regions  of  size  4L0,  with  the  effects 
of  scales  less  than  or  equal  to  2 L0  averaged  out.  This  well  known  block  spin  method 
is  due  to  Kadanoff  (1977). 

The  above  process  can  be  expressed  in  terms  of  a  transformation  T,  which  is 
applied  iteratively:  THq  — >  H\,  TH\  — >  Hi,  TH%  — >  H3 . . ..  At  each  stage,  the  length 
scales  are  changed:  L0  — >  2 L0, 2 L0  — >  4L0  . . .,  and  in  order  to  compensate,  the  spin 
variables  are  also  scaled  in  an  appropriate  fashion  (see  next  subsection).  It  is  this 
rescaling  which  leads  to  renormalization,  and  the  set  of  transformations  {T}  defines 
a  semi-  group:  hence  renormalization  group.  If  iterating  this  transformation  leads  to 
the  result  that  an  integer  N  exists  such  that 

THn  =  Hn+1  =  Hn,  for  all n>N 

then  Hn  =  HN  (say)  is  referred  to  as  a  fixed  point  In  the  case  of  critical  phenomena, 
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this  corresponds  to  a  critical  point.  We  conclude  this  brief  discussion  of  the  magnetic 
case  by  indicating  how  its  methodology  may  be  made  to  come  closer  to  a  continuum 
problem  (as  in  fluid  motion),  which  is  what  concerns  us. 

Suppose  we  make  the  Fourier  transformation  from  x-space  to  wavevector  k-space. 
Further  suppose  that  we  have  eliminated  (in  some  way)  the  smallest  length  scales,  less 
than  or  equal  to  A  1,  say,  so  that  we  are  left  with  a  Hamiltonian  which  is  a  function 
of  continuous  variables  in  k-space,  H( A).  Then  we  may  formally  eliminate  further 
modes  in  a  band  A/6  <  k  <  A,  as  follows  (1  <  b  <  oo).  Denote  the  three-dimensional 
spin  field  by  Sa( k),  where  a  =  1,  2,  or  3,  and  take  the  following  steps: 

1.  Integrate  out  all  Sa(k)  for  which  A/6  <  k  <  A. 

2.  Rescale  the  remaining  modes  of  the  spin  field  by  enlarging  wavevectors  by  a 

factor  6. 

3.  Multiply  each  ,5Q(k)  by  a  constant  factor  (b. 

The  parameter  6  is  known  as  the  spatial  rescaling  factor,  while  (b  is  the  spin  rescaling 
factor,  and  involves  technicalities  of  the  magnetic  case,  which  we  shall  not  pursue 
here.  If  we  write  the  above  three-step  process  as 

H'  =  TbH, 

then  the  set  %,  1  <  6  <  oo  is  called  the  renormalization  group.  It  is  often  remarked 
that,  since  Tb  1  does  not  exist,  this  is  strictly  a  semi-group. 

We  shall  now  consider  RG  as  applied  to  the  2D  Ising  model. 

2.3  Renormalization  group  theory  in  critical  phenomena 

2.3.1  Iterative  Hamiltonian  structure 

As  we  discussed  already,  an  important  and  early  application  of  RG  was  to  critical 
phenomena.  Here  we  shall  provide  a  more  detailed  review  of  the  calculation  of  Wilson 
(1975;  Wilson  and  Kogut,  1974)  for  the  2D  Ising  model.  In  the  2D  Ising  model 
of  ferromagnetism,  one  considers  a  lattice  of  points  and  on  these  lattice  sites  are 
particles  (called  ’spin’  particles)  which  (i)  can  take  on  only  two  values  [  i.e.,  ‘spin  up’ 
or  ‘spin  down’],  and  (ii)  can  only  interact  with  its  nearest  neighbors.  The  system  is  in 
thermal  equilibrium,  with  a  lattice  temperature  For  simplicity,  we  shall  consider 
a  square  lattice  with  many  spin  particles.  ^From  equilibrium  statistical  mechanics 
(Huang,  1963),  it  is  known  that  all  thermodynamic  properties  of  this  system  can  be 
determined  from  the  partition  function 

z  =  Y.  expH<>[s  ;  i  s 
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where  the  Hamiltonian  Ha,  for  the  2D  Ising  model, 

Ho[s]  =  K£X>n.sn+i  (2) 

K  =  —J/kDT,  where  J  is  the  spin-spin  interaction  coupling  constant  [with  J  >  0  for 
ferromagnetism] .  The  sum  in  (1)  is  evaluated  over  all  possible  spin  configurations. 
The  outer  sum  in  (2)  is  over  all  lattice  points,  while  the  inner  sum,  for  the  represen¬ 
tative  lattice  site  vector  n  =  (rq,  n2),  is  a  sum  over  the  four  nearest  neighbors  of  the 
square  lattice 

'y  ]  ■c>n^n+i  ®n^(ni+l,tia)  "b  l,ri2)  b  *^n^(ni,n2+l)  "b  •^ln'^(ni,n2  1)*  (^) 

i 

The  difficuly  in  determining  the  thermodynamic  properties  of  this  microscopic 

Ising  model  is  in  the  evaluation  of  the  summation  in  (1)  :  there  are  just  too  many 
degrees  of  freedom2.  The  important  things  to  note  are  (i)  the  original  problem  in¬ 
volved  only  nearest-neighbor  interactions,  (ii)  the  coupling  strength  K  is  related  to 
the  strength  of  the  interactions,  J,  and  the  lattice  temperature 

2.3.2  Iterative  decimation  of  the  number  of  degrees  of  freedom 

Suppose  that  the  degrees  of  freedom  are  decimated  by  restricting  the  summation  in 
Eq.  (1)  to  run  only  over  half  the  lattice  sites.  Thus  the  partition  function  Z  must 
now  be  determined  from 

Z  =  ^2exp[Hi[(r}}  (4) 

l 

where  the  new  lattice  sites  are  denoted  as  cr,  the  sum  is  over  these  sites  only,  and  Hi 
is  a  Hamiltonian  to  be  determined  so  that  (4)  yields  the  identical  partition  function 
as  Eq.  (1). 

It  can  be  shown  that  the  new  Hamiltonian  Hi  [cr]  takes  the  form 

Hi  [a]  =  (N2/2)A  +  2BJ2  Wk  +  Bj^  fiifijfik  +  C  ^  (5) 

nn  nnn  p 

where  “nn”  denotes  nearest-neighbor  coupling,  “nnn”  next-nearest-neighbor  coupling 
and  ‘p’  the  plaquette  four-spin  interactions  (Hu,  1982)  with 

A{k)  =  ln2  +  (1  /8)  [J,n{cosh4K)  +  4ln(cosh2K)]  (6) 

B(k )  =  (l/8)ln(cosh4K)  :.“j 

C{k)  =  (1/8)  [ln{coshAK)  -  4ln{cosh2K )]  (8) 

2We  are  actually  using  this  2D  Ising  model  as  an  example  of  how  to  treat  problems  that  have 
so  many  degrees  of  freedom  that  new  techniques  must  be  introduced  to  solve  them.  The  2D  Ising 
model  is  actually  exactly  soluble,  but  that  is  not  of  direct  concern  to  us. 
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Here  A  is  simply  a  normalization.  Eq.  (5)  should  be  compared  to  the  original  Ising 
Hamiltonian  H0,  eq.  (2).  First  we  see  the  expected  changes:  (a)  the  nearest  neighbor 
interactions  term  has  strength  2B(K),  rather  than  K,  and  (b)  there  is  a  normaliza¬ 
tion  constant  A(K)  introduced  due  to  the  change  in  the  correlation  length  of  the 
decimated  system,  What  is  of  significant  interest  is  the  appearance  of  new,  non- 
Ising,  interactions:  diagonal  nearest  neighbor  interactions  [strength  B(k)],  and  four 
spin  interactions  [with  strength  C(k)].  Therefore,  although  the  original  Hamiltonian 
has  only  nearest-neighbor  coupling,  other  types  of  interactions  are  generated  by  the 
renormalization  group  to  compensate  for  the  reduction  in  the  number  of  degrees  of 
freedom.  This  is  a  general  feature  of  the  renormalization  group  (Hu,  1982). 

We  now  proceed  to  the  second  iteration  to  determine  the  Hamiltonian  H2  for 
N2 /22  lattice  sites.  H2  is  to  be  determined  so  that  the  partition  function  Z  remains 
invariant  under  decimation  of  spins.  If  all  the  4  terms  are  kept  in  (5)  in  determining 
H2,  then  one  cannot  calculate  the  sum  over  the  decimated  spins  exactly.  To  pro¬ 
ceed,  one  can  start  with  the  nearest-neighbor  interaction  and  treat  the  other  types 
of  interactions  perturbatively.  Wilson  proceeds  under  the  assumption  that  the  4-spin 
interaction  (and  any  higher  order  spin  interaction)  is  weak  and  can  be  neglected,  but 
it  is  critical  that  the  W  interactions  be  retained.  The  self-consistency  of  these 
assumptions  can  be  checked  at  the  end  of  the  calculation. 

Hence,  if  at  the  i  -  Ith  stage,  the  Hamiltonian  is  tfi_i[s],  with  its  dependence 
on  the  Ising  (original)  nearest-neighbor  interaction  strength  Rj_i  and  the  (non-Ising) 
diagonal-nearest-neighbor  interaction  A  i ,  then  at  the  ith  iteration,  the  Hamiltonian 
f/j[s]  is  a  function  of  Ki  and  L; .  It  is  important  to  note  that  only  in  degenerate  cases 
is  the  Hamiltonian  form  invariant  under  RG  transformations. 

2.3.3  RG  transformations 

Now  for  sufficiently  weak  interactions,  one  can  expand  the  interaction  strength  coef¬ 
ficients  in  a  Taylor  series  in  K  =  —J/KbT: 

B(K)  =  ln{cosh4K)/ 8  «  K2  +  0(A4:  :.9j 

C(K)  =  (1/8 ){ln{cosMK)  -  4ln{cosh2K)\  =  0(K 4),  (10) 

allowing  us  to  drop  the  4-spin  interaction  term  in  (5).  Note  further  that  the  coeffi¬ 
cients  for  “uu”  and  “mm”  interactions  are  of  the  same  order  in  K,  requiring  us  to 
keep  these  newly  created  diagonal  “mm”  interactions.  To  determine  the  RG  trans¬ 
formation,  we  must  note  the  following:  (a)  the  diagonal  nearest-neighbor  interaction 
at  the  ith  stage,  Li,  arises  from  the  coupling  coefficient  B(Ki-i)  =  R?_i  +  0(Kf_ i), 
on  \ising  (7):  i.e., 

U  =  KU  (11) 

(b)  the  nearest-neighbor  interaction  at  the  ith  stage  arises  from  2  terms:  (i)  the 
nearest-neighbor  interactions  at  the  i  -  Ith  stage,  with  coupling  coefficient  B(Ki  i)  = 
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Kf_x,  and  (ii)  the  diagonal-nearest-neighbor  interactions  at  the  i  —  1th  stage  which 


yields  a  coupling  coefficient  Lj_i  to  K\.  This  is  because  the  diagonal-  nearest-neighbor 
interaction  at  the  i  -  1th  stage  will  become  the  nearest-  neighbor  interaction  at  the 


jth 


stage  due  to  the  decimation  of  spins. 

Hence  the  RG  transformations  to  proceed  from  the  i  — 

Ith  to  the  ith  stage  are 

Ki  =  2KU  + 

(12) 

h  =  KU 

(13) 

with 

L0  =  0 

(14) 

since  initially  we  are  dealing  with  an  Ising  model  which,  by  definition  only  has  nearest- 
neighbor  interactions. 

The  fixed  points,  Kt  and  L»,  of  (12)-  (14)  are  immediately  determined  from: 


Kt  —  2  K*  +  L*  (15) 

L.  =  K\  (16) 

There  are  3  possible  fixed  points: 

K*i  =  0,  L*i  =  0  (17) 

K„  2  =  oo,  L,  2  =  oo.  (18) 

K*3  =  1/3,  L»3  =  1/9  (19) 


and  the  initial  domain  of  attraction  in  the  K-L  plane  for  these  fixed  points  are  readily 
determined  (with  L0  =  0):  (a)  for  0  <  K0  <  0.39209,  L0  =  0  :  stable  fixed  point 
is  (0,0);  (b)  for  K0  >  0.3921,  L0  =  0  :  stable  fixed  point  is  (oo,  oo);  (c)  for  K$  = 
0.39209,  Lq  =  0  :  unstable  fixed  point  is  (1/3, 1/9) 

In  the  Ising  problem,  the  relevant  fixed  point  is  the  saddle  point  (c).  The  fixed 
point  (a)  corresponds  to  K»  =  0  =  J/kBT, ,  i.e.,  T,  ->  oo.  Physically,  this  corresponds 
to  a  perfectly  disordered  system,  with  only  short  scale  fluctuations  that  tend  to  zero 
as  the  number  of  RG  iterations.  The  fixed  point  (b)  corresponds  to  the  temperature 
T,  =  0,  and  this  corresponds  to  a  perfectly  ordered  system,  again  with  negligible 
fluctuations.  Hence,  as  far  as  critical  phenomena  are  of  interest,  these  fixed  points 
axe  unimportant.  However,  the  fixed  point  (c)  is  important  since  even  after  all  the 
RG  transformations,  there  are  still  many  length  scales  present.  Critical  exponents 
prove  to  be  properties  of  this  fixed  point. 
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3 


Application  of  RG  To  Turbulence 


3. 1  Similarity  and  differences  between  fluids  and  critical  phe¬ 
nomena 

The  macroscopic  motion  of  a  fluid  may  be  thought  of  as  exhibiting  two  ‘phase  tran¬ 
sitions5.  That  is,  as  we  increase  the  Reynolds  number  of  a  given  flow,  we  may  expect 
first  a  transition  to  turbulence  and  second  a  transition  to  scaling  behaviour  of  fully 

developed  turbulence.  In  the  latter  case,  this  means  the  existence  of  an  intermediate 
range  of  wavenumbers  in  which  the  energy  spectrum  takes  the  form  of  a  power-law. 
It  is  this  transition  which  concerns  us  here.  Specifically,  we  will  assume  the  iner¬ 
tial  range  ‘5/3’  law  for  the  energy  spectrum  developed  by  Kolmogorov  [1941].  As 
an  example,  consider  a  stationary  homogeneous  isotropic  turbulence  where  energy 
is  injected  with  a  forcing  spectrum  peaked  a  wavenumber  ko(=  1/Zo),  where  Iq  is 
the  integral  length  scale.  The  energetics  of  eddies  with  wavenumbers  smaller  than 
a  given  wavenumber  k  are  determined  by  an  equilibrium  between  the  injection,  the 
losses  due  to  viscous  dissipation  and  the  energy  flux  function  II(fc)  that  moves  energy 
to  higher  wavenumbers  (i.e.,  to  smaller  scales).  Here,  II  (A;)  depends  upon  the  inertial 
terms  in  the  Navier-Stokes  equation,  and  is  expressible  in  terms  of  the  triple  veloc¬ 
ity  correlation  (see  for  example,  Rose  and  Sulem,  1978).  The  problem  is  defined  by 
three  parameters  plus  the  underlying  Navier-Stokes  equation.  These  parameters  are 
the  kinematic  viscosity,  the  total  rate  of  dissipation,  £',  which  equals  the  energy  flux 
function,  and  the  characteristic  lenth  Z0.  To  make  contact  with  the  critical  phenom¬ 
ena,  the  external  length  scale  l0  might  be  considered  analogous  to  a  lattice  constant 
(Nelkin,  1974). 

It  is  of  course  difficult  to  identify  the  onset  of  scaling  behaviour  as  such  and  a  con¬ 
venient  way  to  think  about  this  problem  is  to  consider  the  case  where  the  Reynolds 
number  (based  on  the  turbulent  microscale,  say)  is  large  enough  for  a  power-law  spec¬ 
trum  to  exist  over  an  appreciable  range  of  wavenumbers.  Now  define  a  local  Reynolds 
number,  using  the  inverse  of  wavenumber  as  the  length  scale.  Then  as  we  scan  through 
the  wavenumbers,  from  the  maximum  towards  smaller  values,  we  are,  in  effect,  in¬ 
creasing  the  local  Reynolds  number  to  the  point  where  power-law  behaviour  begins 
(In  other  words,  the  boundary  between  the  inertial  and  viscous  ranges).  Conceptu¬ 
ally,  therefore,  the  simplest  application  of  RG  to  turbulence  involves  the  progressive 
elimination  of  high-A;  modes,  with  the  fixed  point  corresponding  to  a  renormalized 
viscosity  in  the  inertial  range  (The  equivalent,  in  terms  of  CFD,  would  be  a  subgrid 
effective  viscosity,  in  large-eddy  simulations).  We  shall  carry  out  these  operations  on 
the  Navier-Stokes  equations. 
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3.2  RG  applied  to  the  Navier-Stokes  equations 

We  introduce  the  Fourier  components  of  the  velocity  field  ua( k,  t)  in  the  usual  way, 

uQ(x,t)  =  ^uQ(k,t)elk-x.  (20) 

k 

Then  the  solenoidal  Navier-Stokes  equation  for  incompressible  fluid  motion  takes  the 

form  (e.g.  see  McComb  1990) 

(^  +  v0k2)ua(k,t)  =  XMa01(k)  J d3ju^(j,t)tt7(k  —  j,t),  (21) 

Map-,(k)  =  (2i)-l{k^Dcn{k)  +  k^D^  k)},  (22) 

A*(k)  =  4*-^k  (23) 

A  is  a  book-keeping  parameter,  which  can  be  used  to  facilitate  the  iterative  solution 
of  the  Navier-Stokes  equation  by  expanding  in  powers  of  A.  One  then  sets  A  =  1  at 
the  end  of  the  calculation.  This  procedure  is  well  known  in  statistical  and  many-body 
physics  and  is  usually  referred  to  as  the  A-expansion.  It  is  believed,  however,  that  this 
expansion  is  divergent  for  the  Navier-Stokes  equation  (Orszag  1977)  -  nevertheless  it 
is  utilized  in  all  RG  theories. 

Turbulence  almost  invariably  occurs  because  there  is  a  mean  rate  of  shear  in  a 
fluid.  If  we  wish  to  study  stationary,  isotropic  turbulence,  then  the  artificial  nature  of 
this  problem  requires  us  to  add  hypothetical  stirring  forces  to  the  right-hand  side  of 
the  Navier-Stokes  equations,  in  order  to  maintain  the  turbulence  against  viscous  dis¬ 
sipation.  Denoting  the  random  force  by  fa(k,t),  it  is  usual  to  specify  its  distribution 
as  multivariate  normal  and  to  choose  the  force-force  covariance  to  be 


where 

and 


</o(k,  k,  If))  =  Dap(k)W(k)5(t  -  If),  (24) 

where  W{k)  has  dimensions  of  velocity2 /time,  and  remains  to  be  specified.  Station- 
arity  requires  that  the  rate  at  which  the  stirring  forces  do  work  is  the  same  as  the 
rate  of  viscous  dissipation  e,  or 


4tt  k2W{k)dk 


=  e, 


(25) 


where 


2v0k?  E(k)dk. 


(26) 


It  has  always  been  normal  practice  in  turbulence  theory,  to  choose  W[k)  to  be 
peaked  near  the  origin,  so  that  its  arbitrary  nature  is  only  of  importance  at  low 
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wavenumbers,  and  a  universal  energy  spectrum  can  develop  at  large  wavenumbers. 
More  recently,  the  application  of  RG  to  stirred  fluid  motion  (Forster  et  al,  1977; 
Yakhot  and  Orszag.  1986)  has  introduced  theories  which  depend  strongly  on  the 
choice  of  W(k).  Accordingly,  questions  of  how  W(k)  is  chosen  and  justified  must 
then  be  considered.  In  particular,  in  this  context,  the  question  of  whether  or  not 
W(k)  should  depend  on  some  characteristic  length  scale  (such  as  the  ultra-violet 
cutoff  A)  will  surface  later  as  a  major  issue. 

In  order  to  try  to  carry  out  the  RG  algorithm  given  above  in  Section  2.2,  we 
consider  the  Fourier  components  to  be  defined  on  the  interval  0  <  fc  <  k0)  where  fco 
is  the  largest  wavenumber  present  and  is  of  the  order  of  the  Kolmogorov  dissipation 
wavenumber  (e.g.  see  McComb  1990).  Then,  for  some  ki,  such  that  h  <  k0,  we  filter 
the  velocity  field  at  k  =  ki.  Here,  k^  =  (1  —  r])k0  with  0  <  r}  <  1.  This  may  be 
expressed  in  terms  of  the  unit  step  functions, 


6  (k)  = 

{l 

if  0  <  k  <  ki 
if  ki  <  k  <  k,0 ; 

(27) 

9*(k)  = 

{J 

if  0  <  k  <  ki 
if  k\  <  k  <  k0, 

(28) 

allowing  us  to  define  the  following  useful  filtered  forms: 

ua(k,t)=0  (fc)uQ(k,f), 

«£(M)  =  e+(k)ua(k,t), 

Mafaik)  =  0  (k)Mapy( k), 

Af^(k)  =  e*(k)Ma^{ k). 

Then  substitution  of  these  forms  into  equation  (24)  allows  us  to  decompose  the  Navier- 
Stokes  equation  into  separate  low-fc  and  high-fc  forms,  viz., 


{^  +  v0k2)ua(k,  t ) 

=  AMa/j7(k)  J  d?j  {up (j,  t)uy  (k  —  j,  f) 

+2u0  (j,  t)u+  (k  —  j,  £) 
+u£(j,t)u+(k-j,t)},  0  <  k  <  ki 

(29) 

(^  +  VQk2)u+{k,  t) 

=  A M+pj (k)  J  d?j  {up(3,t)uy  (k-j,  t) 

+2up(j,t)itf(k-  j,t) 

+Up (j, t)u+  (k-j, t)}.  ki  <k<  k0 

(30) 

If  we  now  try  to  carry  out  the  first  of  the  three  steps  outlined  at  the  end  of 

Section  2.2,  for  the  lattice-spin  Hamiltonian,  we  can  adapt  that  procedure  to  the 
present  (Navier-Stokes)  case,  as  follows: 
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1.  Solve  equation  (30)  for  u+. 

2.  Substitute  this  solution  into  equation  (29)  for  the  u  and  do  a  partial  average 
over  the  u+. 

3.  Terms  resulting  from  this  procedure  which  are  linear  in  u~  can  be  interpreted 
as  contributing  an  increment  to  the  turbulent  viscosity. 

The  main  problems  which  arise  due  to  mode-mode  coupling,  are  easily  understood: 

•  First,  the  solution  of  (30)  for  u+  contains  terms  in  u~.  When  these  are  sub¬ 
stituted  into  the  right-hand  side  of  (29),  the  result  is  a  triple  nonlinearity  in 
the  u  ,  which  can  then  generate  even  higher-order  nonlinearities  in  subsequent 
iterations. 

•  Second,  averaging  out  the  high-wavenumber  modes  requires  the  hypothesis 

(u  u+u+)  =  u  (u+u+), 

which  cannot  be  strictly  true,  as  u  and  u+  are  just  parts  of  the  same  velocity 
field  and  are  not  statistically  independent.  <  ...  >  represents  the  partial  average 
over  the  high-wavenumber  modes. 

It  cannot  be  emphasized  too  strongly  that  these  two  problems  are  fundamental 
stumbling  blocks  in  applying  RG  to  the  Navier-Stokes  equation.  It  is  the  way  in  which 
they  are  tackled  which  distinguishes  one  RG  theory  of  turbulence  from  another.  We 
will  discuss  in  detail  two  methods  for  dealing  with  these  problems:  the  recursive 
mode  elimination  and  the  conditional  averaging  method.  In  the  approach  to  RG  of 
Forster  et  al  (1977),  Fournier  and  Frisch  (1993),  which  culminated  in  the  theory  of 
Yakhot  and  Orszag  (1986),  suitable  assumptions  are  made  so  that  these  problems  are 
avoided.  This  approach  is  described  in  the  Sec.  6  of  this  review. 


3.3  Iterative  conditional  averaging:  a  summary 

We  now  introduce  a  method  of  eliminating  turbulent  modes  which  is  based  on  the 
use  of  a  conditional  average  to  distinguish  between  amplitude  and  phase  correlation 
effects.  It  has  its  roots  in  the  method  of  iterative  averaging,  which  was  developed 
over  a  number  of  years  as  a  possible  method  of  applying  the  renormalization  group 
approach  to  real  fluid  turbulence  (McComb  1982, 1986,  1990).  However,  an  essential 
feature  of  the  more  recent  work  is  the  formal  treatment  of  the  conditional  average  and 
the  development  of  methods  of  approximating  its  relationship  to  the  usual  ensemble 
average  (McComb  and  Watt  1990,  1992;  McComb,  Roberts  and  Watt  1992). 

The  basic  idea  is  that  the  turbulent  velocity  field  in  wavenumber  space  may  be 
decomposed  into  two  distinct  fields.  One  is  a  purely  chaotic  field ;  while  the  other  is  a 
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correction  field,  and  carries  all  the  phase  information.  Application  of  this  decomposi¬ 
tion  to  a  thin  shell  of  wavenumbers  in  the  dissipation  range  allows  the  elimination  of 
modes  in  that  shell;  with  the  usual  mode-coupling  problems  being  circumvented  by  the 
use  of  a  conditional  average.  The  (conditional)  mean  effect  of  the  eliminated  modes 
appears  as  an  increment  to  the  viscosity,  with  terms  of  order  rf  being  neglected, 
where  tj  is  a  dimensionless  measure  of  bandwidth  thickness,  such  that  0  <  r)  <  1. 
An  iteration  (with  appropriate  rescaling)  to  successively  lower  shells,  reaches  a  fixed 
point,  corresponding  to  a  renormalized  turbulent  viscosity. 

3.4  Recursive  renormalization  group:  a  summary 

If  we  return  to  the  equilibrium  Ising  model,  in  Sec.  2.3,  in  which  the  Hamiltonian  Ho 
involves  only  nearest-neighbor  interactions,  it  was  shown  by  Wilson  (1975)  that  the 
effective  Hamiltonians  Hn  (n  >  1)  now  involve  both  nearest-neighbor  and  diagonal 
nearest-neighbor  interactions.  Form  invariance  of  Hn  is  then  imposed  by  assuming 
that  the  interactions  of  still  higher  order  are  small.  Similarly,  in  our  form  of  recursive 
RG  (Zhou  et  al.,  1988,  1989;  Zhou  and  Vahala,  1992,  1993a, b),  as  applied  to  the 
Navier-Stokes  equations,  form  invariance  can  be  imposed,  when  both  the  quadratic 
and  the  RG-induced  cubic  nonlinearities  are  retained,  by  truncating  the  A-expansion. 

The  recursive  RG  theory  extends  Rose’s  (1977)  treatment  of  the  linear  problem  of 
passive  scalar  convection  to  Navier-Stokes  turbulence  and  eliminates  the  small  scales 
recursively.  Navier-Stokes  equation  include  (1)  a  renormalized  eddy  viscosity  and  (2) 
a  triple  product  in  the  fluid  velocity.  The  eddy  viscosity  is  calculated  by  means  of  a 
difference  equation  (Zhou  et  al.  ,  1988,  1989)  or,  in  the  limit  77  — ►  0,  by  means  of 
a  differential  equation  (Zhou  and  Vahala,  1992).  The  resulting  eddy  viscosity  term 
exhibits  a  mild  cusp  behavior  in  the  renormalized  momentum  equation  (Zhou  et  al., 
1988,  1989).  The  triple  velocity  products  are  absent  only  when  there  is  a  spectral 
gap  between  the  subgrid  and  resolvable  scales  (Zhou  et  al.,  1988).  Beyond  the  second 
iteration,  the  RG-induced  triple  products  also  contribute  to  the  eddy  viscosity  in  the 
renormalized  Navier-Stokes  equation.  The  recursivee  RG,  therefore,  has  potential  for 
capturing  a  variety  of  dynamical  features  that  depends  on  the  interplay  between  local 
and  nonlocal  triad  interactions  (Zhou,  1993a, b;  Zhou  et  al.  ,  1995). 

It  is  of  particular  interest  to  determine  the  effect  of  the  new  triple  velocity  product 
on  the  resolvable  scale  energy  transfer  and  the  difference,  if  any,  from  that  of  the 
usual  Navier-Stokes  quadratic  velocity  product,  and  we  do  this  by  considering  their 
individual  contributions  to  the  eddy  viscosity.  Within  recursive  RG  theory  (Zhou  and 
Vahala,  1993a),  it  has  been  shown  that  the  triple  velocity  product  in  the  renormalized 
momentum  equation,  which  produces  a  fourth  order  velocity  product  in  the  energy 
equation,  removes  energy  from  the  resolved  scales.  The  averaged  fourth-order  product 
can  be  decomposed  into  a  product  of  averaged  second-order  products  and  formally 
takes  the  form  of  a  gradient  diffusion  process  with  eddy  viscosity  vT.  The  triple 
term  contribution  to  the  eddy  viscosity  is  zero  as  k  — >  0  and  increases  rapidly  as 
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k  —>  kc.  Here  kc  is  the  cutoff  wavenumber  which  separates  the  resolvable  and  subgrid 
scales.  The  spectral  eddy  viscosity  is  simply  the  sum  of  the  contributions  from  the 
momentum  equation  and  the  triple  nonlinear  term  and  it  appears  to  be  in  qualitative 
agreement  with  the  closure  theory  (Kraichnan,  1976;  Chollet  and  Lesieur,  1981)  and 
direct  numerical  measurement  (Domaradzki  et  al.,  1987;  Lesieur  and  Rogallo,  1989; 
Zhou  and  Vahala,  1993a).  In  particular,  it  predicts  the  correct  asymptotic  behaviors 
of  the  eddy  viscosity  as  k  — ►  0  and  k  — >  kc.  The  methodology  has  been  applied  to 
the  passive  scalar  being  advected  by  incompressible  turbulence  (Zhou  and  Vahala, 
1993b). 
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Part  III 

Technical  Aspects  of  selected  RG 
approaches 

4  Iterative  conditional  averaging:  details 

In  this  section  we  amplify  our  earlier  summary  of  the  method  and  show  in  some  more 
detail  how  the  iterative  conditional  averaging  technique  is  used  the  implement  the 
general  strategy  of  the  renormalization  group  for  Navier-Stokes  fluid  turbulence. 

4.1  The  conditional  average 

First  of  all,  we  formulate  the  operation  of  taking  a  conditional  average.  Then  we  use 
this  average  in  such  a  way  that  we  can  eliminate  the  triple  nonlinearity  referred  to 
earlier. 

The  idea  is  quite  simple.  Putting  it  at  its  most  basic  level,  we  select  from  the  full 
ensemble  of  turbulence  realizations  a  sub-ensemble,  the  members  of  which  have  their 
low- A;  modes  equal  to  ua  (k,  t).  Then  we  perform  averages  of  functions  of  the  (k,  t ) 
over  this  sub-ensemble. 

However,  there  is  more  to  it  than  this.  If  u(k,  t)  is  the  solution  of  the  Navier- 
Stokes  equation,  corresponding  to  prescribed  boundary  conditions,  then  we  are  faced 
with  (in  principle)  a  deterministic  process  and  to  prescribe  u-  is  to  prescribe  u+. 
That  is  to  say,  if  u“  is  invariant  under  our  conditional  average,  then  so  also  is  u+.  In 
order  to  get  round  this  problem,  we  invoke  the  defining  characteristic  of  deterministic 
chaos.  This  is  to  the  effect  that  any  uncertainty  in  the  specification  of  the  system 
will  be  amplified  exponentially;  so  that  as  time  goes  on  the  difference  between  almost 
identical  solutions  will  increase  to  the  point  of  unpredictability.  In  the  present  case, 
we  replace  the  concept  of  time  going  on,  by  the  number  of  steps  of  the  cascade 
in  wavenumber.  That  is,  our  ideas  about  the  turbulent  cascade,  and  particularly 
ideas  about  localness  of  energy  transfer,  suggest  that,  if  we  prescribe  conditions  at 
wavenumber  Ay,  then  u+(k0,f)  will  be  unaffected,  provided  that  A:0  is  much  larger 
than  k\.  In  other  words,  as  the  bandwidth  becomes  large  (j]  —>  1),  the  conditional 
average  of  u+  becomes  free  of  constraint  and  we  can  expect  that 

(u+(k0,  t))c  -►  (u+(k0,t)). 

On  the  other  hand,  for  k0  — ►  ki,  it  is  intuitively  clear  that  the  conditional  average 

must  tend  to  become  effectively  deterministic,  with 

(it+(k0,  t))c  u+(ko,t). 
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We  shall  return  to  this  at  a  later  stage,  but  for  the  present  it  should  be  appreciated 
that  the  point  at  issue  here  is  one  of  phase  correlation.  If  two  realizations  are  slightly 
out  of  phase  at  wavenumber  k\ ,  then  we  can  (given  the  nature  of  turbulence)  expect 
that  this  phase  difference  will  amplify  exponentially,  so  that  phase  correlation  will 
decline  throughout  the  band  as  one  moves  from  ki  to  k0.  Hence,  for  w+(k0,  t)  to  be 
chaotic,  despite  the  prescription  n+(k1,t)  =  u  (ka,t),  it  is  a  requirement  that  the 
bandwidth  be  large  enough.  This  imposes  a  lower  bound  on  possible  values  of  the 
bandwidth  parameter  Tj- 

We  may  take  account  of  these  various  aspects  by  introducing  a  fuzzy  criterion 
for  our  conditional  average.  We  now  choose  as  our  sub-ensemble,  the  subset  of  real¬ 
izations  for  which  the  low- A:  modes  differ  from  u“(k ,t)  by  a  small  amount  $“(k,  t). 
Obviously  $  is  an  arbitrary  criterion  (apart  from  the  constraint  that  u  +  $  must 
be  a  possible  solution  of  the  Navier-Stokes  equation)  and  should  be  chosen  to  satisfy 
the  conditions  which  we  wish  to  impose  upon  our  average. 

We  shall  not  go  into  this  procedure  any  further  here,  and  in  order  to  simplify  the 
algebra,  we  shall  omit  the  4  from  the  equations.  The  interested  reader  who  wishes 
to  pursue  the  matter  will  find  a  full  account  in  the  papers  by  McComb  and  Watt 
(1992)  and  McComb  et  al  (1992). 

4.2  Conditionally-averaged  equations  for  high  and  low  wavenum¬ 
bers 

Let  us  now  denote  the  operation  of  taking  a  conditional  average  over  the  modes  in 
the  band  k\  <  k  <  k,0  by  angle  brackets,  with  a  subscript  0.  This  notation  per¬ 
mits  the  subsequent  generalization  to  subscripts  1,2 , ,n,  as  we  remove  shells  of 

wavenumbers  progressively.  Then,  we  list  the  ideal  defining  properties  of  the  condi¬ 
tional  average  as 

(u  (k,  t))o  =  u  (k,  t),  (31) 

(u  (k,t)u  (k',t))0  =  u  (k,t)u  (k ',*),  (32) 

and  so  on,  for  products  of  the  low-wavenumber  modes  of  any  order. 

We  now  conditionally  average  both  equations  (29)  and  (30)  which  are,  respectively, 
the  low- A;  and  high-A;  filtered  Navier-Stokes  equations. 

First,  we  obtain  the  conditional-averaged  NSE  on  the  interval  0  <  k  <  ki,  by 
averaging  according  to  equations  (31)  and  (32),  to  obtain 

(£  +  ^oA;2)w“(k ,t)  =  J  d3j  {u^  (j,  t)u~ (k  —  j,  t) 

+2u^  (j,  t)(u+  (k  —  j,  t))a 

+(n^(j,t)w+(k- j,t))0}.  (33) 
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We  now  repeat  the  steps  just  taken,  but  this  time  we  apply  them  to  equation  (33) 
for  the  high-fc  modes.  Thus,  we  get 

( £  +  vok2)u+( k,  t)  =  M+^(k)  j  d3j  {up  (j,  (k  -  j,  t) 

+2 (j,  t)u+  (k  -  j,  t)  +  uj (j,  (k  -  j,  0},  (34) 

and  taking  the  conditional  average  of  each  term  according  to  equations  (32)  and  (33), 
gives 


+i/0fc2)(«a(M))o  =  M+p^(k)  J d3j  {up(j,t)u^(k  — j,t) 

+2up  (j ,  t)  (u+  (k  —  j ,  f )  )o 
+(u|(j,t)u+(k- j,t))0}. 

(35) 

4.3  First-shell  elimination  using  the  two-field  decomposition 

Our  objective  now  is  to  solve  equation  (35)  for  the  conditional  average  (u| (j,  t)u+  (k  —  j,  t))0 
and  substitute  the  result  back  into  (33),  in  order  to  have  a  closed  equation  for  the  low- 
wavenumber  modes.  In  order  to  do  this,  we  shall  ultimately  have  to  reckon  with  the 
need  to  relate  conditional  averages  to  full  ensemble  averages.  Accordingly,  we  begin 

this  section  with  the  two-field  decomposition  which  is  our  basis  for  this  procedure. 

Let  us  write  the  exact  decomposition: 

u£  (k,  t)  =  v+ (k,  t)  +  A+  (k,  t) ,  (36) 

where  o+(k,  t)  is  any  other  realisation  of  our  turbulent  ensemble.  In  other  words, 
u+(k,  t)  has  exactly  the  same  statistical  properties  as  «£(k,  t),  but  has  no  phase 
relationship  to  u~(k,  t).  It  follows,  by  definition,  that  A+(k,  t)  is  simply  a  measure 
of  the  phase  difference  (in  the  band  of  modes  to  be  eliminated)  between  the  two 
realizations.  It  also  follows  that,  from  the  point  of  view  of  the  realization  under 
study,  u+(k,  t)  is  the  purely  chaotic  part  of  the  field  and  A+(k,  t)  is  the  correction 
field  which  carries  all  the  phase  information. 

We  are  now  in  a  position  to  write  down  an  expression  relating  the  conditional 
average  of  the  high-wavenumber  part  of  the  velocity  field  to  its  full  ensemble  aver¬ 
age.  Taking  the  conditional  average  of  both  sides  of  equation  (33),  it  may  be  shown 
McComb  et  al  (1992),  or  it  is  intuitively  obvious,  that 

(ut  (k>  <))o  =  (k,  t))  +  (A+  (k,  t))0.  (37) 

Now  we  seek  a  relationship  between  v+  and  u+,  which  is  such  that  the  conditionally 
averaged  correction  term  (A"*")o  may  be  neglected  as  small.  In  other  words,  we  need 
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an  ansatzior  the  correction  term  and  naturally  this  will  depend  on  the  physical  nature 
of  the  system  that  we  are  studying. 

In  the  case  of  macroscopic  fluid  turbulence,  we  are  guided  by  the  well  established 

idea  that  turbulent  energy  transfer  in  wavenumber  takes  the  form  of  a  cascade  and  is 
therefore  to  some  extent  local  in  wavenumber.  In  terms  of  our  present  approach,  we 
take  this  to  mean  that,  in  any  particular  realisation,  the  effect  on  phase  correlation 
of  the  mode-mode  coupling  is  short-range  in  k- space.  Thus,  on  a  statistical  picture, 
based  on  many  such  realisations,  modes  which  are  widely  separated  may  be  taken  to 
be  independent  of  each  other.  Hence,  providing  that  the  bandwidth  parameter  77  is 
not  too  small,  we  can  assume  that  u+  (k0 ,  t)  is  independent  of  u+(ki,t),  in  the  sense 
that  we  can  write 

(u+(k0,  <))o  =  (u+(k0,  *))  =  (v+(ko,  t)),  (38) 

where  the  last  step  follows  from  the  definition  of  v+,  as  another  realization  of  the 
turbulence  ensemble  with  the  same  statistical  properties  as  u+,  but  with  no  phase 
relationship  to  u  . 

This  now  leads  us  towards  a  natural  ansatz  for  the  relationship  between  v+  and 
u+.  Relying  on  the  fact  that  we  are  dealing  with  a  problem  in  continuum  mechanics, 
we  take  (k,  t)  to  be  given  by  a  first  order  truncation  of  the  expansion  of  u+  (k,  t ) 
in  Taylor  series  about  k  =  k0,  thus: 

i£(M)  =  Wo(k0,t)  +  (k  —  ko).Vfctia(k,  t)  |fc=fco  +0(tj2).  (39) 

Note  that  we  conclude  that  terms  of  order  rj1  have  been  neglected  because  the  maxi¬ 
mum  value  of  |k  —  ko|  is  Tjk^.  Hence  it  follows  that  we  have 

(A+(k,  t))o  =  0(rf).  (40) 


It  also  follows,  therefore,  that  we  are  simultaneously  imposing  both  Tipper  and 
lower  bounds  on  acceptable  values  of  77.  On  the  one  hand,  77  must  be  large  enough  for 
us  to  assume  that  u+(k0,t)  is  independent  of  u+(k1,t);  while,  on  the  other  hand,  77 
must  be  small  enough  for  us  to  neglect  terms  which  are  of  order  77 2  in  equation  (39). 

Then,  with  all  these  points  in  mind,  above  equations  yield  for  the  viscosity  acting 
on  the  explicit  scales: 

Vi  =  Vo  +  5v0,  (41) 


where  the  formula  for  the  increment  to  viscosity  is 

=  ±  +0(rT), 

k?  J  u0j 2  +  i/0|k-j|2 


(42) 


with  0  <  k  <  ki,  ky  <  j,  |k  —  j|  <  fcQ  and  Q+  is  merely  an  extension  of  the  spectral 
density  to  the  v+  field.  The  coefficient  L(k,  j)  is  given  by 


L(k,j)  =  — 2M^(k)M/3,«(j)£>i7(|k  — j|) 

_  [jti(fc2  +  j2)  -  kj{  1  +  2ju2)](l  -  fj?)kj 
k2  +  j2  —  2k, j n 


(43) 
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where  fj,  is  the  cosine  of  the  angle  between  the  vectors  k  and  j. 

We  extend  the  procedure  to  further  wavenumber  shells,  as  follows: 

(a)  Set  ua  (k,  t)  =  ua(k,  t)  in  equation  (19),  so  that  we  now  have  a  new  NSE  with 
effective  viscosity  V\(k)  for  Fourier  modes  on  the  interval  0  <  k,  <  ki. 

(b)  Make  the  decomposition  of  (61),  but  this  time  at  k  =  fc2,  such  that  u+(k,  t)  is 
now  defined  in  the  band  <  k  <  k\. 

(c)  Repeat  the  procedures  used  to  eliminate  the  first  shell  of  modes  in  order  now 
to  eliminate  modes  in  the  band  k^  <  k  <  k\. 

In  this  way,  we  can  progressively  eliminate  the  effect  of  high  wavenumbers  in  a 

series  of  bands  kn+\  <  k  <  kn ,  where 


K  =  (l-v)nko  5  o  <  77  <  1,  (44) 

with,  by  induction,  the  recursion  relation  for  the  effective  viscosity  given  by 

Vn+l{k)  =  vn{k)  +  Sun(k ),  (45) 


where  the  increment  of  order  n  takes  the  form 

r  (U\  -  1  f  ^3  •Mk)j){(3(0|j=fcn  +  ~  fcn)^p|i=fc„  +  0(r]2)} 

n[  )  ~  k?J  dJ  »n(j)j2  +  M |k— j|)|k  — j|2 

Also,  we  may  form  an  energy  equation  for  the  explicit  scales,  hence  obtaining  the 
renormalized  dissipation  relation,  viz: 


[kn  2 un{k)E{k)dk  =  e,  (47) 

Jo 

which  may  be  compared  with  the  unrenormalized  where  kn  and  vn  should  be  replaced 
by  k,Q  and  u0,  respectively. 

If  we  now  assume  that  the  energy  spectrum  in  the  band  is  given  by  a  power  law 
and  make  the  scaling  transformation 


kn+i  hkn, 


(48) 


where,  for  compactness,  we  define  h  =  (1  -  j?),  it  follows  from  equations  (45)  and  (46) 
that  the  effective  viscosity  may  be  written 

v{k,nk!)  =  a1/2e1/3k~4/3vn{k!)  (49) 


where  a  is  the  constant  of  proportionality  in  the  assumed  spectrum.  Now  the  recur¬ 
sion  relation  becomes 


with 


P„+i (k')  =  hi/3un(hk')  +  h  4/Hi>n(k!) 


5vn(k/)  = 


i  f  t3  7  H KTW 

47r/c/2  J  J  vn(hj')j'2  +  vn(hl')la 


(50) 

(51) 
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(52) 


for  the  wavenumber  bands  0  <  k!  <  1;  1  <  j',l'  <  h  1  where  1/  =  k;  —  j;  | ,  and 

Q’  =  /j11/3  —  — /i14/3(f  —  hr1)  +  h.o.t. 

3 

Iteration  of  equations  (50)  and  (51)  reaches  a  fixed  point  and  once  this  fixed  point 
is  found  we  can  calculate  the  Kolmogorov  constant  by  solving  equations  (47)  and 
(49)  simultaneously.  One  merit  of  taking  a  as  a  test  is  that  it  does  have  known 
experimental  values,  albeit  scattered  in  the  range  1.2  <  a  <  2.2.  Our  calculated  value 
of  the  Kolmogorov  spectral  constant  is  a  =  1.60  ±  0.01  independent  of  bandwidth  in 
the  range  0.25  <f]<  0.45,  and  in  good  agreement  with  experiment.  For  values  of  rj 
outside  this  range,  the  calculated  a  diverges  from  the  experimental  value.  At  large 
values  of  77,  this  is  due  to  the  breakdown  of  the  first-order  Taylor  series  approximation, 
while  at  small  values,  one  is  seeing  the  effects  of  mode  coupling,  which  would  invalidate 
the  assumption  that  u(k0)  is  independent  of  u(ki). 


5  Recursive  RG 

5.1  Introduction 

We  shall  apply  recursive  renormalization  group  (RG)  procedures  to  the  problem  of 
subgrid  modeling.  Subgrid  modeling  is  necessary  for  the  high-Reynolds  number  tur¬ 
bulent  flows  of  interest  because  of  the  limitations  of  current  and  forseeable  super¬ 
computers.  A  motivation  is  that  the  spectral  transport  coefficients  (such  as  the  eddy 
viscosity)  determined  from  recursive  RG  theory  can  be  compared  to  those  arising 
from  closure-based  theories  (Kraichnan,  1976;  Chollet  and  Lesieur,  1981).  It  should 
be  noted  that  the  transport  coefficients  in  these  closure  theories  are  determined  over 
the  entire  resolvable  scales,  and  are  a  function  of  k  in  the  resolvable  scales. 

In  particular,  we  point  out  here  that  in  e-RG,  a  small  parameter  e  is  introduced 
through  the  forcing  correlation  function.  Yakhot  &  Orszag  (1986)  have  to  extrapo¬ 
late  from  e  <C  1  to  e  — >  4  in  order  to  reproduce  the  Kolmogorov  energy  spectrum. 
Furthermore,  it  is  also  necessary  to  take  the  distant  interaction  limit,  k  — >  0.  Thus,  it 
is  difficult  to  compare  the  wave-number  dependent  transport  coefficients  (Kraichnan, 
1976;  Leslie  and  Quarini,  1979;  Chollet  and  Lesieur,  1981),  with  that  determined 
from  e-RG. 

In  recursive  RG,  no  attempt  is  made  to  introduce  a  special  form  of  overlapping 
as  in  conditional  averaging,  but  one  proceeds  directly  with  standard  averaging  and 
handles  the  triple  nonlinearity  directly.  The  basic  differences  between  the  recursive 
and  e  RG  procedures  are  that  in  recursive  RG: 

(i)  The  e-expansion  is  not  applied. 

(ii)  The  turbulent  transport  coefficients  are  determined  for  the  whole  resolvable 

wavenumber  scales, 
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(iii)  Higher  order  nonlinearities  are  generated  in  the  renormalized  momentum 

equation  and  play  a  critical  role  in  determining  the  transport  coefficients. 

(iv)  RG  rescaling,  as  in  conditional  averaged  RG,  is  performed. 

It  should  be  emphasized  that  there  are  two  singular  limits:  h  — ►  1  and  k  — *•  0.  A 
careful  analysis  must  be  done  regarding  these  two  limits  and  the  associated  averaging 
operations.  We  will  address  this  issue  here  in  the  present  review. 

5.2  Renormalized  momentum  equation 

For  notation  consistency,  we  shall  use  “+,  instead  of  “>,  <”  which  are  typically 
employed  in  recursive  RG. 

The  first  and  third  terms  on  the  RHS  of  (29)  are  symmetric  in  j  and  |k  —  j |  in 
terms  of  their  respective  wavenumber  constraints  in  wavenumbers.  As  a  result,  the 
distant  interaction  limit  k  — ►  0  has  no  effect  on  the  existence  of  these  terms,  and 
these  terms  will  give  rise  to  the  standard  quadratic  nonlinearity  (first  term  of  Eq.  29) 
and  eddy  viscosity  (third  term  of  Eq.  29).  However,  the  second  term  on  the  RHS  of 
(29)  has  the  following  constraint:  j  is  in  the  subgrid  while  |k  —  j |  is  in  the  resolvable 
scales.  Specifically,  the  consistency  condition  requires  that,  for  small  k,  j  satisfies 

j  >  ki  and  j  <  ki  +  kz 

where  k  -  j  =  kjz.  Since  \z\  <  1,  the  range  of  integration  must  be  0(k). 

Thus,  the  second  term  on  the  RHS  of  Eq.  (29)  can  not  contribute  in  the  limit 
k  — *  0  since  the  integrand  is  bounded.  Now  it  is  well  known  that  the  higher  order 
nonlinearities  axe  induced  by  this  second  term  under  discussion.  Since  this  term  is 
absent  in  the  k  — *■  0  limit,  we  conclude  that  the  higher  order  nonlinearies  will  not 
contribute  to  the  renormalized  momentum  equations  and  recursion  relation  for  the 
transport  coefficients  in  the  distant  interaction  limit,  k  — >  0.  However,  they  will 
contribute  to  the  renormalized  Navier-Stokes  for  0  <  k  <  k\.  Ater  obtaining  the  final 
renormalized  Navier-Stokes  equation 


[d/dt  +  v(k)k?]ua( k,  t )  =  /Q(k,  t)  +  Ma/j7(k)  J  d3  jup( j,  t)uy( k  -  j,  t) 

+2 - ^-HrMa/g7(k)  f  d3jd3j'-^jMp^ 8V  (j)u^(j',  t)uY{j  -  j',  t)u7(k  -  j,  t).  (53) 

v(kc)kc 3  J  J  3 

by  removing  other  subgrid  shells  iteratively,  one  can  write  down  the  corresponding 
wavenumber  restrictions  and  perform  the  same  analysis  on  the  k  — >  0  limit. 

5.3  Galilean  invariance  of  the  renormalized  Navier-Stokes 

Here,  we  txirn  our  attention  to  the  question  of  the  Galilean  invariance  of  the  renormal¬ 
ized  Navier-Stokes  equations.  The  importance  of  Galilean  invariance  in  turbulence 
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modelling  has  been  emphasized  by  Speziale  (1985).  To  be  consistent  with  the  basic 
physics,  it  is  required  that  the  description  of  the  turbulence  be  the  same  in  all  inertial 
frames  of  reference.  The  appearance  of  the  triple  nonlinear  term,  which  is  a  function 
of  the  resolvable  scales  velocity  fields,  makes  the  property  of  the  Galilean  invariance 
of  our  recursive  RG  procedure  not  transparent.  We  now  show  that  the  renormalized 
Navier  -Stokes  equation  is  Galilean  invariant  (Zhou  and  Vahala,  1993b). 

The  Galilean  transformation  is 


x  — >  x*  —  U0<*  t  -►  t*. 

Here  U0  is  an  uniform  velocity  field.  Thus,  one  has 


u 


u*-U0, 


—  =  —  d  _  d  d 

dx  dx*  dt  dt*  0  dx* 


While  the  Galilean  transformation  for  the  Navier-Stokes  equation  in  physical  space 
is  trivial,  the  Galilean  transformation  in  wavenumber  space  is  less  obvious,  due  to 
the  lack  of  differential  operations.  For  convenience,  we  first  review  how  Galilean 
invariance  is  preserved  for  the  Navier-Stokes  equation  in  the  wavenumber  space. 

Under  the  Galilean  transformation,  the  LHS  of  the  Navier-Stokes  equation  [cf. 
Eq.  (21)]  becomes 


+  UopikpU*a(k* ,  t)  +  v0k*2[-U0a5(k*)  +  n*  (k*,  t)] 

=  +  Uopi^uX k*,  t)  +  v0k*2u*a(k*,  t ) 

where  in  the  last  step,  we  have  used  the  the  6  function  property  k*2S(k*)  —  0. 

Also,  under  the  Galilean  transformation,  the  RHS  of  the  Navier-Stokes  equation 
becomes 


MQ/3l(k*)  / d3;?K(j*,t)  -  UQf}6(n}[u}(k*  -  j *,i)  -  EM(k-  -  j*)] 

=  MQf3l{k*)  J d3j*u}(j\t)u;(k* -i\t)  +  iU^u*a{k\t) 
where  we  have  used  the  property  of  the  5  function,  the  incompressible  condition,  and 

Mah(k*)UO0U;(k*,t)  =  Ut D0  kpU*a  (k* ,  f)  /  (2i) 

Thus,  as  expected,  the  Navier-Stokes  equation  is  invariant  under  a  Galilean  transfor¬ 
mation  due  to  the  cancellation  of  the  second  term  on  the  RHS. 

To  show  that  the  renormalized  Navier-Stokes  equation  is  invariant  under  a  Galilean 
transformation,  we  need  only  consider  the  recursive  RG  induced  triple  nonlinear  term, 
denoted  by  N St' 


24 


NSt  =  2Mafh(k)  J  d3 jd3 j'(~)i/3 Q  -  j'l  t)uY(j',  t)u7(k  -  j,  t)  (54) 

It  is  important  to  note  that  j  is  in  the  subgrid. 

Under  a  Galilean  transformation,  Eq.  (54)  becomes 


nsi = -j">] 

[«;.()“,  t)  -  -  j*,  t)  -  rh„sik-  -  j-)] 

Since  j*  is  in  the  subgrid  scale,  while  f*  and  k*  are  in  the  supergrid,  d(k*  -  j*) 
and  5(j*  —  j'*)  can  never  be  simultaneously  satisfied.  As  a  result, 


(55) 

Now  only  one  term  in  Eq.  (55)  could  violate  the  Galilean  invariance  of  the  renormal¬ 
ized  Navier-Stokes  equation.  However, 

J  d3r5(j'*)u},(r -r,t)  -  t£(jv) 

This  is  not  permissable  since  u*j,  =  u—*p  and  j*  is  restricted  to  the  subgrid.  Thus 
NSt  =  NSt-  Hence  the  triple  term  is  Galilean  invariant. 


NS*  =  2Ma^(k*)  j d3j*d3f*(Z-)V3- 


(j 


(k*)j* 


5.4  The  effect  of  cubic  nonlinearity 

We  consider  the  contribution  of  the  triple  nonlinear  term  in  the  renormalized  mo¬ 
mentum  equation  to  the  eddy  viscosity.  The  second  moment  for  the  velocity  field  is 
defined  as 

Uafi(. k,  t)  =<  ua( k,  *)up(-k,  t)  >  .  (56) 

The  time  evolution  of  Uap( k,  t)  is 

dUaj^  =  -2v(k)k*Ual3(k,  t)  +  7S(  k,  t)  +  7j(  k,  t).  (57) 

In  this  equation,  T^(k,  t)  is  the  standard  energy  transfer  from  the  quadratic  nonlin¬ 
earity.  In  contrast,  TjQ(k,t)  =  -2vT{k)k?  E{k)  is  the  energy  transfer  arising  from  the 
RG  induced  triple  nonlinearity.  It  is  readily  shown  that 
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vr{k)  = 


(58) 


Tja(k,t)  =  1  1  fk+kcd-d  ^(fc»J»g)lk~j|  "  y+1/3 

2E(k)k?  -  2v(kc)  k?  Jkc  3  Z  i/(k-j) 

It  has  been  shown  that  that  vT(k)  is  the  major  contributor  to  the  cusp-like  be¬ 
havior  of  the  spectral  eddy  viscosity  as  k  — >  0  (Zhou  and  Vahala,  1993a).  This  type 
of  term  was  found  to  be  the  major  contribution  to  the  strong  cusp  in  the  spectral 
eddy  viscosity  found  from  the  closure  models  (Kraichnan,  1976;  Leslie  and  Quarini, 
1979;  Chollet  and  Lesieur,  1981). 

Rose  (1977)  discussed  the  role  of  the  triple  nonlinear  terms  in  physical  space.  He 
pointed  out  that  it  represents  the  possibility  of  an  exchange  of  scalar  eddies  between 
the  resolvable  and  subgrid  scales.  This  effect  is  an  inherent  property  of  measurements 
made  on  the  passive  scalar  system  with  instruments  which  have  a  spatial  resolution 
limited  to  an  eddy  width  size  greater  than  1  / kc. 

5.5  Difference  equations  for  the  renormalized  eddy  viscosity 

After  the  removal  of  the  (n  +  l)th  subgrid  shell,  the  spectral  eddy  viscosity  in  the 
renormalized  momentum  equation  is  determined  by  the  recursion  relation 


where 


and 


vn+ 1(*0  =  vn{k)  +  Si/n{k) 


Svn(k)  =  §  £  / 

h= 0 


.  £(fc,j,g)lk-j|  v 
V(j)rMlk-  j|)lk-j|2’ 


(59) 

(60) 


L{k,j,q )  =  — 


kj{\  -  z2)[zq2 
Q 2 


(61) 


with  k  •  j  =  kjz  and  q  =  |k  -  j|..  This  difference  equation,  after  rescaling,  has  been 
solved  by  Zhou  et  al.  (1988,  1989)  and  fixed  points  were  readily  determined  for 
finite  h  <  0.7.  Note  that  the  spectral  eddy  viscosity  shows  a  mild  cusp  as  k  -»  kc, 
in  qualitative  agreement  that  that  of  closure  theory.  However,  it  was  very  difficult 
to  determine  fixed  points  for  finer  subgrid  partition  factor  h  >  0.7.  In  the  next 
subsection  we  shall  pass  to  the  differential  subgrid  limit  h  —>  1  and  determin  an 
ordinary  differential  equation  (o.d.e)  for  the  renormalized  eddy  viscosity  over  the 
entire  resolvable  scale  which  can  be  readily  integrated. 


5.6  Differential  equations  for  the  renormalized  eddy  viscosity 

The  differential  limit,  h  ->  1,  is  a  singular  one  and  this  point  has  been  discussed 
recently  (Zhou  and  Vahala,  1992).  In  particular,  it  is  related  to  the  assumption  of 
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local  versus  non-local  interactions  in  k.  In  this  section  we  will  calculate  the  eddy 

viscosity  tinder  the  differential  equation  limit  for  recursive  RG. 

For  recursive  RG  we  will  find  that  the  differential  equations  hold  throughout  the 
resolvable  wavenumber  range  0  <  k  <  kc.  This  should  be  contrasted  with  the  e  -  RG 
eddy  viscosity  differential  equation  which  is  valid  only  in  the  k  — >  0  limit. 

5.6.1  The  differential  equation  limit,  h  — >  1 

We  now  derive  the  differential  equation  from  which  the  transport  coefficients  for  finite 

k,  0  <  k  <  1,  are  determined.  The  o.d.e.  in  the  distant  interaction  limit  (k  =  0) 
will  be  derived  in  the  next  subsection.  After  the  rescaling,  we  rewrite  the  recursion 
relation  in  the  form 


vn+1  (k)  -  h{y+1V3vn(hk)  =  h^3Sun{hk).  (62) 

For  h  — »  1,  the  number  of  interation  n  — >  oo.  Similarity  consideration  leads  to 

vn+i{k)  ->  v{k),  n  ->  oo.  (63) 

The  LHS  of  Eq.  (62)  becomes 

v(k)  -  [1  -,)]  -  „[^  +  +OW]  (64) 

As  noted  earlier,  the  partial  average  of  Rose  must  be  employed  in  order  to  ensure 
the  existence  of  the  differential  limit.  The  partial  average  is  introduced  since  the 
distinction  between  the  resolvable  and  subgrid  scales  becomes  fuzzy  in  the  limit  of  a 
differential  subgrid  partitioning,  h  — »■  1. 

Following  Rose,  we  first  change  the  variable  from  j,  z  to  j,q  ,  with  djdz  = 
(qfkj)djdq,  so  that  the  RHS  of  Eq.  (62)  becomes 


Sv[k)=jimx)— 

VL  -dq 


L(k,j,  q) 


+ 


LiKj,q)  (l)fo+i)/3 


kj'v{j)k?v(  |k-j|)|k 
L(k,\,q) 


•  j]2|k  —  j|^  kc 


h<q<i+k  k3v2(l)qy 


f  ,  2)/3 

VJi<j<i+k  3  kW(  1)  3 


where 

L(k,  1,  q)  =  k(  1  -  z2)(k  -  zq*)/q2 

and 


L(k,j,  1)  =  kj(  1  -  z2)(kj  -  z). 


(65) 


(66) 

(67) 
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As  a  result,  the  fixed  point  renormalized  eddy  viscosity  v(k)  is  determined  from 
the  o.d.e.  at  0(A) 

k^dlr + !Lr,'(k)  =  w>[Mk) + B-(k)]  (68) 

where 


dq 


qy  i 


(69) 


BJk)  =  1  f*1  djL(k,  j,  l)/»  2>p  (70) 

Here,  z  is  evaluated  at  j  =  1  and  q  =  1,  respectively  in  the  L(k,j,  q)  expression.  The 
o.d.e.  for  the  momentum  equation  eddy  viscosity  is  readily  solved  (Zhou  and  Vahala, 
1993a). 


5.6.2  Differential  equations  in  the  k  —*  0  limit 

In  the  limit  k  — »  0,  we  have  seen  that  the  triple  nonlinearities  induced  by  RG  do  not 
contribute  to  the  eddy  viscosity.  As  a  result,  the  recursion  relation  will  now  contain 
only  the  usual  quadratic  contribution.  We  further  simplify  the  analysis  by  taking 
the  standard  subgrid  linear  propagator  G^(|k  -  j|)  =  +  W)»(|k  —  j|)j  ~  Gr^1(|j|)  as 

k  — >  0. 

The  limits  of  the  integration  are  given  by 


1  <  j  -  kz  <  1//,  1  +  kz  <  j  <  1  //  +  kz 

(71) 

Thus,  one  has 

5vn{k)  =  S  - 

-F-G 

(72) 

where  the  integral  limits  for  these  terms 

are 

/•i//  r1 

y,  *y,* 

for 

S 

(73) 

/»!  r\-\-kz 

L  dzL  d} 

for 

F 

(74) 

f0  pl/f 

/  cfe  /  dj 

J-l  Jl/f+kz 

for 

G. 

(75) 

Terms  F  and  G  are  the  corrections  to  the  symmetric  term  S.  They  are  important 
for  a  finite  bandwidth  /.  However,  it  is  easy  to  show  that  F  +  G  =  0  for  /  — *  1  in 
the  k  — >  0  limit.  Hence 


m  bsiL  /', dz[1  -  "2| Iz+ -y  - ')]  =  4^ 
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while  the  LHS  of  the  same  equation  yields 


km+y±lm^v±l 


v(k), 


as 


A;  — >  0 


since  is  bounded  as  k 
Thus,  as  k,  — >  0, 


0. 


v(k  — ♦  0) 


3  8  1 

y  +  115  i/2(l)' 


(77) 


Again,  the  o.d.e.  for  the  momentum  equation  eddy  viscosity  is  readily  solved.  We 
observe  that  the  eddy  viscosity  has  a  similar  plateau  structure  as  k  — »  0.  As  k  — »  kc, 
eddy  viscosity  displays  a  weak  cusp  like  behavior  as  k  — >  k,c.  In  this  case,  those  curves 
are  similar  to  that  of  Zhou  et  al  (1988;  1989). 


5.7  Recursive  RG  Spectral  eddy  viscosity 

The  spectral  eddy  viscosity  is  simply  the  sum  of  the  contributions  from  the  momentum 
equation  and  that  of  the  effect  of  the  RG  induced  triple  nonlinear  term  in  the  energy 
equation.  It  is  apparent  that  our  calculation  is  in  qualitative  agreement  with  that 
from  the  closure  theory  (Kraichnan,  1976;  Chollet  and  Lesieur,  1981).  and  direct 
numerical  measurements  (Damaradski  et  al.,  1987;  Lesieur  and  Rogallo,  1991;  Zhou 
and  Vahala,  1993).  In  particular,  it  predicts  the  correct  asymptotic  behaviors  of  the 
eddy  viscosity  as  k  — >  0  and  k  — >  kc  (Kraichnan,  1976). 

5.8  Numerically  evaluated  eddy  viscosity 


These  conclusions  can  be  tested  directly  using  numerical  simulation  databases.  In¬ 
deed,  energy  transfer  and  spectral  eddy  viscosity  can  be  analyzed  using  results  from 
direct  numerical  simulations  by  introducing  an  artificial  cut  at  a  wavenumber  kc  that 
is  smaller  than  the  maximum  resolved  wavenumber  km  of  the  simulation.  With  this 
fictitious  separation  between  the  subgrid  and  resolvable  scales,  it  is  possible  to  eval¬ 
uate  the  effect  of  the  subgrid  kc  <  k  <  km  on  the  resolved  scales  k  <  kc.  We  form  an 
energy  equation  from  the  momentum  equation  and  introduce  the  following  notation: 
T><(k)  and  T>:>(k)  represent  the  spectrum  of  energy  transfer  to  mode  k  resulting 
from  interactions  with  one  and  both  modes  above  the  cutoff  kc  respectively.  Measure¬ 
ments  of  numerical  simulation  databases  indicate  the  following  (Zhou  and  Vahala, 
1993a): 

•  T>>(k)  removes  energy  throughout  the  resolvable  scales  in  a  manner  consistent 
with  the  notion  of  eddy  viscosity. 

•  T><(k.)  removes  energy  from  the  last  resolved  octave  that  was  transferred  there 

by  the  resolved  scale  transfer;  that  is,  it  allows  the  local  flow  of  energy  through  kc. 
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It  is  the  most  important  subgrid  effect  near  kc  and  accounts  for  most  of  the  energy 
flow  from  the  resolved  scales. 

The  subgrid  spectral  eddy  viscosity  v>>(k)  and  v><(k)  can  be  determined  from 
T>>{k)  and  T>K(k)  for  a  given  energy  spectrum,  E(k).  Specifically,  u>>(k)  = 
—T>>(k)/2k2E(k)  and  v><{k)  =  —T><(k)/2k2E{k).  Two  important  features  of 
the  quadratic  contribution  v>>{k)  should  be  stressed.  First,  its  positive  constant 
asymptote  at  small  k  indicates  that  the  concept  of  modeling  this  contribution  as  an 
eddy  viscosity  in  analogy  to  the  molecular  viscosity  is  plausible,  and  second,  its  value 
decreases  monotonically  as  k  increases  toward  kc.  This  indicates  that  if  we  include 
only  the  contribution  of  quadratic  velocity  products,  there  is  no  eddy  viscosity  cusp  at 
the  cutoff  kc.  The  most  important  feature  of  v>K(k)  is  the  sharp  increase  at  k  — *  kc. 

6  Reconsideration  of  the  YO  theory 

Robert  Rubinstein 

Institute  for  Computer  Applications  in  Science  and  Engineering 
NASA  Langley  Research  Center 
Hampton,  VA  23681 


The  first  comprehensive  application  of  RG  methods  to  turbulence  was  the  Yakhot- 
Orszag  (1986)  theory  (YO).  This  theory  built  upon  earlier  work  by  Forster  et  al 
(1977)  and  Fournier  and  Frisch  (1983),  but  completed  the  development  suggested  by 
its  predecessors  with  impressive  successes  inchiding  simple  analytical  calculations  of 
the  Kolmogorov  constant  3  and  several  constants  of  interest  in  turbulence  modeling 
(Speziale,  1991).  Nevertheless,  this  theory  continues  to  generate  controversy.  Three 
points  will  be  addressed  in  this  section  about  which  particularly  strong  objections 
have  been  raised: 

1.  YO’s  treatment  of  mode  elimination 

2.  the  “correspondence  principle” 

3.  the  e-expansion  and  distant  interaction  limit 

This  section  draws  heavily  on  work  of  Woodruff  (1992,  1993)  which  emphasizes  the 

connections  between  the  YO  theory  and  Kraichnan’s  (1959)  direct  interaction  ap¬ 
proximation  (DIA).  Although  this  viewpoint  tends  to  de-emphasize  the  importance 
of  mode  elimination  characteristic  of  renormalization  group  theories,  it  is  consistent 
with  the  remark  of  Eyink  (1994)  that  YO  is  not  a  true  RG  theory  in  any  case. 

In  this  Section,  the  following  notation  will  be  used: 

k  =  (k,  fi)  k  =  (p,  u)  q=  { q,  H  —  u>) 

3The  Kolmogorov  constant  had  been  computed  earlier  by  Kraichnan  (1964)  by  numerical  inte¬ 
gration  of  the  LHDIA  closure. 
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where  the  triangle  condition 


k  =  p  +  q 


holds  throughout. 

6.1  YO’s  treatment  of  mode  elimination 

It  was  noted  at  the  end  of  Sect.  3.2  that  iterated  mode  elimination  leads  to  two 
fundamental  problems:  it  generates  nonlinearities  of  higher  order  than  appear  in 
the  Navier-Stokes  equations,  and  it  requires  some  simplifying  assumptions  to  permit 
averaging  over  the  high  wavenumber  components  of  the  motion.  These  difficulties 
arise  because,  by  itself,  mode  elimination,  like  the  field  theoretic  functional  formula¬ 
tions  of  the  turbulence  problem,  does  nothing  more  than  reformulate  the  equations 
of  motion.  In  order  to  make  progress,  statistical  hypotheses  must  be  introduced. 
This  circumstance  is  not  surprising:  for  example,  equilibrium  statistical  mechanics 
requires  Gibbs’  hypothesis,  a  very  strong  assumption  which  certainly  does  not  follow 
from  Newtonian  mechanics. 

Nevertheless,  subsequent  work  has  shown  that  the  original  treatment  of  mode 
elimination  by  YO  was  not  entirely  satisfactory.  For  example,  their  proposal  that 
higher  order  nonlinearities  are  negligible  in  a  certain  perturbative  sense,  was  contra¬ 
dicted  by  Eyink  (1994).  4  But  by  placing  the  YO  theory  in  the  setting  of  DIA-like 
closures,  these  problems  are  understood  by  invoking  the  statistical  hypotheses  of  these 
closures.  Thus,  higher  order  nonlinearities  can  indeed  be  considered,  but  in  the  con¬ 
text  of  higher  order  versions  of  DIA  (Martin  et  al,  1978).  Accordingly,  the  absence 
of  these  nonlinearities  in  the  YO  theory  merely  reflects  the  order  of  approximation 
chosen,  and  requires  no  further  justification. 

The  treatment  of  averaging  over  the  small  scale  motion,  which  appears  to  proceed 
as  if  motions  of  different  scale  were  independent,  has  been  discussed  in  Kraichnan’s 
(1959)  original  presentation  of  DIA:  the  same  requirement  arises  in  the  derivation  of 
the  DIA  response  equation.  Heuristically,  the  DIA  closure  assumes  that  the  velocity 
field  is  only  weakly  non-Gaussian.  Therefore,  the  motions  of  different  scales,  which 
are  uncorrelated  because  of  the  kinematic  hypothesis  of  homogeneity,  are  independent 
to  leading  order.  This  statistical  hypothesis  also  justifies  the  breakup  of  fourth  order 
moments  into  products  of  second  order  moments.  A  theory  in  which  the  motions 
of  different  scale  are  not  independent  to  leading  order,  but  in  which  fourth  order 
moments  are  treated  by  the  quasi-Gaussian  hypothesis,  must  be  carefully  formulated 
to  avoid  inconsistency. 


4YO  had  asserted  that  higher  order  nonlinearities  are  ‘irrelevant’  as  this  is  understood  in  Wilson’s 
(1974)  theory  when  e  =  0,  although  they  are  only  marginal  when  e  =  4.  Eyink  (1994)  demonstrated 
that,  on  the  contrary,  higher  nonlinearities  are  marginal  regardless  of  t. 
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6.2  The  correspondence  principle 


The  starting  point  of  YO’s  analysis  is  the  Navier-Stokes  equations  driven  by  a  random 
force: 

-inui(k)  -  Mimn (k)  /  dpdq  um(p)un(q)  =  fo(k)  (78) 

J k=p+Q 

where  the  Gaussian  random  force  /j  is  characterized  by  its  correlation  function 


<  /,(£)/#')  >=  2D(2it)d+'kr*Dij(k)8(k  +  k!) 


(79) 


This  force  is  white  noise  in  time.  The  exponent  y  is  treated  as  a  variable  for  purposes 
of  the  subsequent  e-expansion,  in  which  e  =  4  +  y  —  d  and  d  =  3  is  the  dimension  of 
space.  The  analysis  leads  to  a  Kolmogorov  spectrum  when  y  =  3  or  e  =  4;  this  is 
therefore  the  case  of  physical  interest.  YO’s  conclusion  is  that  the  nonlinear  term  the 

Navier-Stokes  equation  is  replaced,  in  the  limit  of  infinite  Reynolds  number,  by  the 

combination  of  random  forcing  by  ft  and  a  scale  dependent  viscosity  v(k),  so  that 

-iQui(k)  +  v(k)k2Ui(k)  =  fi(k )  (80) 


Mode  elimination  is  used  to  obtain  the  recurrence  relation 

<fr(*)  ■  D 

dk  v2k? 


(81) 


where  the  constant  A  is  computed  from  the  theory. 

Whereas  it  is  generally  agreed  that  Eq.  (78)  provides  a  plausible  model  of  isotropic 
turbulence  provided  the  random  force  /  is  concentrated  at  large  scales  and  therefore 
provides  an  energy  source,  the  introduction  of  a  force  acting  on  all  inertial  range 
scales  appears  to  lack  fundamental  justification.  YO’s  model  Eqs.  (79)  and  (80)  can 
be  compared  to  the  DIA  Langevin  model  (Kraichnan,  1976) 

-iQu(k)  + 1 j(k)Ui(k)  =  fi(k )  (82) 


where  the  damping  function  and  force  correlation  are  expressed  in  terms  of  the  DIA 
response  and  correlation  functions  G  and  Q  by 

rj(k)  =  2iMrmn(k)  [  dpdq  Dmr (p) Dns (q) G(p) Q (q)  (83) 

Jk=p+q 

F(k)=<fi(ic)fi(k')  >  /S(k+k')  =  —AMirnn (k) Mjrs (k)  f  dpdq  Dmr(p)Dns(q)Q(p)Q(q) 

J  k=p-\-q 

(8A) 

Eq.  (82)  is  a  generic  model  in  statistical  mechanics  which  replaces  the  effects  of  an 
infinity  of  nonlinear  interactions  on  any  one  mode  by  a  random  force  acting  against 
a  generalized  damping;  DIA  applies  this  description  to  a  problem  which  is  far  from 
thermal  equlibrium.  The  “fixed  point”  RG  model  Eqs.  (79)  and  (80)  formally  resem- 
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bles  the  DIA  Langevin  equation  model.  But  the  damping  v(k)k?  in  in  the  RG  model 
Eq.  (79)  is  Markovian,  so  that 

V(k)=v(k)  (85) 

only  and  the  forcing  in  Eq.  (80)  is  white  noise  in  time,  so  that 

F(k)  =  F(k )  (86) 

only.  Neither  of  these  conditions  holds  for  the  DIA  Langevin  model. 

To  investigate  the  connection  between  these  models,  write  following  Kraichnan 

Q(p)  =  Q(p)R(p )  (87) 

where  R  is  the  time  correlation  function.  Perform  the  frequency  integration  in  Eq. 
(84)  and  evaluate  the  result  in  the  long  time  limit  in  which  =  0.  This  limit 
corresponds  to  observing  the  system  over  times  long  compared  to  any  characteristic 
correlation  time  of  the  true  DIA  random  force.  The  result  is 

F(k)  =  -4Mimn(k)Mjrs(k)  [  dpdq  Dmr(p)Dns(q)Q(p)Q{q)Q(k,p,  q)  (88) 

Jk=p+q 


where 


&(k,p,q )  =  f  duR(p)R(q)  |n=o 

J  OO 


In  this  limit,  the  random  force  is  white  in  time.  Further,  in  Kolmogorov  scaling, 


Q(Xp)  =  X  n/3Q(p) 
&(Xk,Xp,Xq)  =  A  2/3Q(k,p,q) 


(90) 

(91) 


consequently,  the  scaling  dimension  of  the  random  force  is  found  to  be  -3: 

F(Xk)  =  X  ~3k,  (92) 


Formally,  in  the  long  time  limit,  the  random  force  in  the  DIA  Langevin  model  has 
the  same  space-time  correlation  as  the  force  postulated  at  the  outset  by  YO. 

It  should  be  noted  that  the  power  counting  which  leads  to  Eq.  (92)  is  purely 
formal,  since  the  actual  force  correlation  integral  in  Eq.  (84),  like  the  integral  of  Eq. 
(83),  is  infrared  divergent  when  evaluated  for  an  infinite  Kolmogorov  inertial  range. 
We  recall  that  these  divergences  actually  cancel  in  the  DIA  energy  equation;  however, 
the  assumption  of  a  -3  force  in  DIA  requires  a  priori  infrared  regularization. 

That  the  -3  force  is  natural  in  the  context  of  any  steady  state  far  from  equilibrium 
with  a  constant  flux  of  some  invisicid  invariant  is  also  suggested  by  the  derivation 
(Rubinstein,  1994a)  of  Bolgiano  scaling  inertial  range  for  buoyant  turbulence  by  ap¬ 
plying  the  YO  formalism  with  forcing  of  the  temperature  equation  only.  From  this 
point  of  view,  the  introduction  by  Lam  (1992)  of  a  distinguished  infrared  scale  in  the 
RG  force  is  debatable:  it  corresponds  to  a  loss  of  the  locality  of  the  inertial  range 
postulated  by  Kolmogorov. 


6.3  The  e-expansion  and  the  distant  interaction  approxima¬ 
tion 

The  e-expansion  is  the  subject  of  an  especially  large  number  of  re-evaluations  and  re¬ 
considerations,  among  which  are  Ronis  (1987),  Lam  (1992)  and  Wang  and  Wu  (1993). 
In  YO’s  original  presentation,  the  e-expansion  is  an  expansion  about  a  logarithmically 
divergent  theory.  An  interesting  alternative  was  suggested  by  Carati  (1990a),  who 
suggested  expanding  about  a  theory  with  vanishing  energy  transfer  (Fournier  and 
Frisch,  1978).  Here,  this  expansion  will  be  considered,  following  Woodruff  (1992),  as 
an  approximation  in  DIA. 

To  complete  the  transition  from  the  DIA  Langevin  model  to  the  YO  theory,  further 
approximations  are  required.  They  are 

•  (a)  evaluate  the  DIA  integrals  in  the  distant  interaction  limit  in  which  k/p,k/q—> 
0 

•  (b)  Markovianize  the  damping 

•  (c)  introduce  an  infrared  cutoff  so  that  the  integrals  in  Eqs.  (83), (84)  are 
restricted  to  p  >  k  and  q>  k  only 

•  (d)  evaluate  the  amplitudes  using  the  e-expansion 

It  has  been  emphasized  by  Woodruff  that  these  approximations  axe  closely  related. 

First,  as  noted  by  Kraichnan  (1987),  the  e-expansion  is  an  expansion  about  a  theory 
in  which  distant  interactions  are  dominant;  accepting  this  point  provisionally,  we 
outline  how  the  distant  interaction  limit  brings  about  the  Markovianization  of  the 
damping  and  forcing. 

Let  the  real  function  T/(£),  0  <  £  <  oo  satisfy 

poo 

H{ 0)  =  1,  H(£)  <  1  for  £  >  0,  /  H(£)d£  <  oo  (93) 

J  0 

Then  standard  properties  of  delta  functions  imply 

A H(X(t  -  s))  ~  5{t  -  s)  for  A  — >  oo  (94) 

Rewrite  Eq.  (83)  in  the  time  domain,  and  evaluate  the  wavevector  integrals  in  the 
distant  interaction  approximation  in  which  k  — *  0,  p,  q  — »  oo.  Then 

f]( k,  t,  s )  =  [  dpdq  B( k,  p,  q )G(p,  t,  s)Q(q,  t,  s ) 

«'k=p-|-q 

~  J dp{km-^-(k, p,p )G(p,t,s)Q(p,t,s) 
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dQ 

-B( k,  p,  p )G(p,  t,  s)kmpmp  x-^(p,  «)} 

where  B(k,  p,  q)  denotes  the  product  of  projection  operators  in  Eq.  (83).  Assuming 
time  stationary  similarity  forms 

G(p,t,s)  =  G(pr(t  —  s)) 

Q(p,t,s)  =  R(pr(t  —  s))Q(p) 

the  properties  Eqs.  (93)  of  H  may  reasonably  be  postulated  of  the  product  GR. 
Therefore,  Eq.  (94)  implies  that  in  this  limit  the  damping  is  Markovian 

i]{k,t,s)  =  6{t-s)i](k) 

and  the  DIA  response  equation  implies  that  the  Green’s  function  is  exponential, 

G(k,  t,  s )  =  exp  [(s  —  t)»7(k)]  for  t  >  s  (95) 

Likewise  evaluating  the  force  correlation  Eq.  (84)  in  the  distant  interaction  limit 
implies  that  the  forcing  is  white  noise  in  time: 

<  /«(k,  *)/,(: k',  a)  >=  5(t-  s)6( k  +  k')f^(k)  (96) 

Computing  the  two-time  correlation  function  from  the  relation 

Qij(k,  t,  .s)<S(k  +  k') 

=  [  dn  G(k,t,r1)  [  dr2  G(k',t,r2)x  <  > 

J  0  JO 

using  Eqs.  (95),  (96)  shows  that  the  fluctuation  dissipation  relation 

Q(k,t,s)  =  Q(k)[G{k,t,s)  +  G(k,8,t)]  (97) 

expressing  the  time  dependence  of  the  correlation  functions  in  terms  of  the  response 
function  is  also  valid  in  this  limit. 

These  simplifications  of  DIA  permit  analytical  evaluation  of  the  inertial  range  con¬ 
stants.  Although  values  of  these  constants  could  be  inferred  from  numerical  solutions 
of  DIA,  say  for  decaying  turbulence,  it  is  natural  to  attempt  analytical  evahiation  as 
well. 

Introducing  Eqs.  (95)  and  (96)  with  the  Kolmogorov  scaling  forms 

E{p)  =  CK£V3k-5/a  (98) 

V(P)  =  CDe1/3p2/3 

into  the  DIA  response  equation  integrated  over  all  time  separations, 

C2  r  P  U//3 

=  +  dpdq  2iMrmn  (k) Dmr (p)  Dns (q)  ^2/3  _|_  ^273^2/3  (^9) 
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Integrating  the  single  time  equation  for  the  correlation  function  with  respect  to 
wavenumber  k  leads  to  a  second  equation,  which  for  Kolmogorov  scaling  gives 

S-  =  -1904  (100) 

This  method  of  evaluating  the  inertial  range  constants  CD  and  CK  fails  because  of 
the  well-known  divergence  of  Eq.  (99)  at  low  wavenumbers.  The  simplest  infrared 
regularization  which  is  consistent  with  Kolmogorov  scaling  is  to  restrict  the  region 
of  integration  to  triads  satisfying  p  >  ak.  Values  of  Cd{ol)  and  Ck{ol)  have  been 
tabulated  by  Leslie  (1972). 

The  e-expansion  can  be  considered  as  a  method  of  infrared  regularization  by  an¬ 
alytic  continuation.  Namely,  continue  to  assume  Eq.  (97),  but  replace  Eq.  (98)  by 
the  general  form 

E(p)  =  CKD2/3k}-2tf3  (101) 

The  scale  independence  of  the  integrated  response  equation  demands 

T](p)  =  CDD^3k2  £/3  (102) 

The  units  of  D,  consistent  with  Eq.  (79),  make  these  equations  dimensionally  correct. 
Substituting  these  scalings  in  the  integrated  response  equation  gives  the  e-dependent 
form  of  Eq.  (99), 

C'2  r  n-1-2<E/3 

=  Jk_p+(i  hpdq  2 iMrmn (k) Dmr  (p )Dns (q)  ^  ^  7/3yj^  e/3  ( ^ 93) 

The  integral  in  Eq.  (103)  is  ultraviolet  divergent  when  e  <  0  and  is  logarithmic  when 
e  =  0.  Woodruff  observes  that  it  is  reasonable  to  evaluate  Eq.  (103)  for  e  >  0  by 
asymptotic  expansion  about  e  =  0.  This  expansion  greatly  simplifies  the  integration 
since  the  ultraviolet  divergence  for  e  <  0  implies  that  the  integral  is  dominated  by 
distant  interactions,  namely  by  wavevector  triangles  such  that  p,  q  — >  oo.  In  this 
limit,  a  simple  analytical  evaluation  of  the  integrals  is  possible.  The  calculation  gives 

C}  =  lA(e)  =  A  :+7lo  +  A1e  +  ---  (104) 

Ck  e  e 

where 


The  constant  A_i  is  distinguished  since  it  is  the  only  one  in  the  series  Eq.  (104)  which 
has  been  evaluated  exactly  in  two  senses.  First,  increasing  the  number  of  “loops,” 
that  is,  considering  terms  in  the  perturbative  solution  of  the  Navier-Stokes  equations 
with  a  larger  number  of  force  correlations,  will  correct  Ap  only  for  p>  0.  It  can  also 
be  shown  (Rubinstein,  1994b)  that  even  at  the  one  loop  level,  correcting  the  distant 
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interaction  approximation  by  power  series  expansions  in  k/p  also  only  corrects  Ap  for 
p  >  0.  Accordingly,  it  is  reasonable  to  evaluate  Eq.  (104)  by  taking  the  leading  term 
only.  Setting  e  =  4, 

Cl=3_ 

CK  20 

which  is  easily  shown  to  be  equivalent  to  YO’s  calculation. 

It  is  sometimes  claimed  that  YO  evaluate  the  inertial  range  constants  by  setting 
e  to  four  and  zero  at  different  places  in  the  same  equation.  Therefore,  it  must  be 
emphasized  that  in  this  calculation,  e  is  never  set  to  any  value  but  4 ■  The  analytical 
procedure  which  leads  to  Eq.  (104)  is  entirely  routine:  it  is  the  evaluation  of  the 
leading  term  in  an  asymptotic  expansion,  not  a  novel  procedure  unique  to  YO. 

The  e-expansion  was  described  earlier  as  an  infrared  regularization  necessary  to 
evaluate  the  right  side  of  Eq.  (99),  which  diverges  when  p  — *•  0.  Triads  with  p  ~  0,  q  ~ 
k  correspond  to  sweeping  of  modes  with  wavevector  |  k  =  k  by  modes  of  much  larger 
scale.  The  dynamic  significance  of  this  divergence  has  been  elucidated  by  Kraichnan 
(1982).  This  divergence  is  removed  in  YO,  and  indeed  in  all  renormalization  group 
approaches  by  focusing  exclusively  on  interactions  for  which  p,q>  k. 

In  fact,  the  e  expansion  is  constructed  so  that  when  e  =  0,  the  dominant  inter¬ 
actions  actually  are  the  distant  interactions  for  which  p,  q  — >  oo:  in  this  case,  the 
integral  in  Eq.  (99)  is  logarithmically  divergent  in  this  limit.  However,  as  Woodruff 
(1993)  notes,  the  integral  becomes  infrared  divergent  when  e  =  3,  and  the  analytic 
continuation  from  e  =  0  to  e  =  4  in  the  YO  theory  becomes  problematic.  Thus,  al¬ 
though  it  is  satisfying  to  be  able  to  compute  the  inertial  range  constants,  and  even  to 
obtain  satisfactory  values  by  a  straightforward  computation,  the  fact  remains  that  the 
analytic  continuation  which  underlies  the  calculation  requires  justification.  Moreover, 
Woodruff  also  suggests  that  one  might  attempt  an  e  expansion  about  this  infrared 
divergence.  Not  unexpectedly,  the  results  are  quantitatively  unsatisfactory,  but  this 
possibility  suggests  that  the  expansion  about  e  =  0  is  not  the  only  one  possible. 

Another  objection  to  this  procedure  can  be  raised  in  connection  with  Eq.  (100): 
the  constant  has  been  obtained  by  exact  evaluation  of  the  triangle  integrals  making 
neither  the  e-expansion  nor  the  distant  interaction  approximation.  However,  the 
integral  can  be  shown  to  be  ultraviolet  divergent  for  e  <  4  and  logarithmic  exactly 
when  e  =  4.  Thus,  there  is  no  possibility  of  an  e  expansion  for  this  integral,  which 
must  be  evaluated  exactly. 
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7  Appendix:  Brief  Description  of  RG  Publications 
Relevant  to  YO 

7.1  The  work  of  Kraichnan 

Kraichnan  (1987)  provides  an  extensive  discussion  on  YO.  First,  he  pointed  out  the 
logical  distinction  between  two  procedures:  the  multiple  scale  elimination  and  the 
e-expansion.  First,  the  evaluation  of  eddy  viscosity  by  perturbative  elimination  of 
successive  small  spherical  shells  of  high-wavenumber  modes  can  be  carried  out  by 
specifying  an  actual  energy  spectrum  (Rose,  1977)  rather  than  introducing  a  forcing 
spectrum.  Kraichnan  discussed  next  the  method  of  “e-expansion”.  In  this  procedure, 
the  properties  of  a  E(k)  ~  k1  2e/3  spectrum  range  are  examined  through  an  expansion, 
in  powers  of  e,  about  the  properties  of  a  spectrum  E{k)  ~  k. 

Kraichnan’s  analysis  is  based  on  estimation  of  the  qualitative  nature  of  eddy 
damping  in  a  spectrum  of  the  form 

E{k)  oc  k3  {kc  <k,<  kd ).  (105) 

The  corresponding  eddy  damping  is  given  by  RG  as 

v{k\p)  oc  {p/k)~^3  (0  <  e  <  4),  (106) 

where  s  =  1  —  e.  For  the  Kolmogorov  inertial  range  spectrum,  e  =  4,  the  eddy 
damping  is  local  for  any  positive  e  since  v[k\p)  — >  0  as  p  —>  oo.  In  YO,  the  explicit 
calculation  of  eddy  damping  is  first  made  in  the  near  neighborhood  of  this  reference 
spectrum.  Now  distant  interactions  really  are  weakly  dominant  for  E(k)  ~  k.  At 
s  —  1  or  e  =  0,  Eq.  (106)  is  replaced  by 

v(k,\p)  oc  ln(kc/p).  (107) 

The  results  are  then  mapped  to  the  Kolmogorov  spectrum  by  taking  e  — >  4.  In 
the  meantime,  YO  retained  terms  through  first  order  in  e,  whereas  4  is  not  a  small 
number.  Hence,  the  convergence  properties  of  the  expansion  are  unclear. 

Kraichnan  (1987)  also  stressed  that  this  type  of  estimation  can  be  done  completely 
by  dimensional  analysis  (Fourier  and  Frisch,  1978;  1983;  Kraichnan,  1982).  The  only 
element  of  the  Navier-Stokes  equation  involved  in  these  discussions  are  the  overall 
energy  conservation  by  nonlinear  terms  and  the  coefficients  of  interaction  of  distant 
wave  vector  triads. 

7.2  The  work  of  Teodorovich  and  Wang  8z  Wu 

Teodorovich  (1987, 93, 94)  used  the  field-theoretical  method  to  re-evaluate  the  results 
of  YO.  We  will  discuss  two  issues  discussed  in  his  most  recent  publication  in  1994. 
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The  first  issue  investigated  by  Teodorovich  is  a  reassessment  of  the  eddy  viscosity 
calculation  of  YO.  This  is  exactly  the  same  subject  examined  in  Wang  k  Wu  (1993). 
Recall  that  in  YO  the  parameter  e  is  used  as  an  small  pertuterbation  expansion  and 
is  an  important  parameter  in  YO’s  eddy  viscosity  expression.  YO  then  took  the 
limit  e  — >  0  when  they  compute  the  coefficient  but  set  e  =  4  in  the  powerlaw  expo¬ 
nent.  However,  both  Teodorovich  and  Wang  k  Wu  (among  many  others  workers) 
noted  that  there  is  no  mathematical  justification  for  doing  this,  and  the  choice  of 
e  =  4  leads  to  unacceptable  results.  Furthermore,  they  found  an  algebraic  error  in 
YO.  Specifically,  in  the  integration  over  wavevector  q,  YO  introduced  a  substitution 
q  — >  k  +  q/2.  However,  YO  did  not  take  this  substitution  into  account  in  the  cor¬ 
responding  transformation  in  the  associated  domain  of  integration.  Correcting  this 
algebraic  error,  they  found  an  eddy  viscosity  which  is  independent  of  the  parameter 
e.  Hence,  the  eddy  viscosity  in  Teodorovich  and  Wang  k  Wu  is  independent  of  e, 
but  has  the  numerical  value  of  the  eddy  viscosity  as  YO  who  took  the  e  — *  0  limit. 
Without  e  dependence  in  the  eddy  viscosity,  Teodorovich  and  Wang  k  Wu  propose 
that  a  good  agreement  with  experimental  data  can  be  obtained. 

It  is  important  to  note  the  limit  e  ->  0  in  YO  was  not  only  used  to  obtain  an 
acceptable  value  for  the  eddy  viscosity,  but  the  limit  e  — >  0  also  plays  an  important  role 
in  eliminating  higher  order  nonlinearities  in  the  RG  procedure.  However,  Teodorovich 
and  Wang  k  Wu  do  not  address  the  closure  problem  if  e  =  4  is  maintained  throughout 
the  analysis. 

The  discussion  now  leads  naturally  to  the  question  on  how  the  local  and  nonlocal 
interactions  are  being  treated  by  YO  -  a  second  major  issue  in  the  work  of  Teodor¬ 
ovich.  For  a  forcing  correlation  which  is  a  function  of  e,  an  energy  spectrum  can  be 
deduced  with  an  e  dependence.  As  discussed  already,  Kraichnan  (1987)  found  that 
the  interactions  are  local  when  e  — ♦  0  but  nonlocal  when  e  =  4.  Hence  the  results 
computed  at  one  limit  may  not  be  used  at  another  limit.  Teodorovich  claims  that 
e  — »  0  corresponds  to  the  dominance  of  local  interaction.  An  analytic  continuation 
in  e  from  a  pole-type  singularity  at  e  =  0  to  the  point  e  =  4  means  neglecting  the 
nonlocal  interactions  with  the  large-scale  modes  (Teodorovich,  1994).  His  analysis  is 
not  transparent. 

7.3  The  work  by  Ronis 

Ronis  (1987)  analyzed  a  model  for  randomly  stirred  fluids  using  RG  on  a  path-integral 

representation  of  the  Navier-Stokes  equation.  Unlike  YO,  he  found  that  a  choice  of 
the  random  force  correlation  exponent  (y  =  -1.5851  in  three  dimensions)  is  needed 
to  give  the  Kolmogorov  5/3  law  at  high  wavenumber. 

Of  particular  interest,  Yuan  and  Ronis  (1992)  discussed  the  issue  on  the  generation 
of  random  force.  They  noted  that  YO  and  Ronis  (1987)  ignored  the  actual  generation 
of  the  turbulence,  e.g.,  at  the  boundaries  of  the  system,  and  the  precise  nature  of  the 
stirring  force  is  not  clear.  In  particular,  they  found  that  there  was  no  a  priori  theory 
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for  the  exponents  used  to  characterize  the  random-force  autocorrelation  function. 
First  of  all,  for  incompressible  turbulence, 

^ t]  +  u(r,  t)  ■  Vu(r,  t)  -  V[p/p)  +  i/0V2u(r,  t)  =  F(r,  t)  (108) 

where  the  force  F  results  from  the  interactions  with  the  boundary  and  as  such  is 
not  stochastic  in  nature.  They  took  the  view  that  a  statistical  force  arises  from  the 
random  force  which  represents  the  effective  force  felt  at  smaller  length  scales  that 
results  from  turbulent,  but  deterministic,  motion  at  larger  scales  as  ‘transmitted’  by 
the  nonlinear  terms  in  NSE.  They  introduced  a  projection  operator  V  such  that 

u'(r,  t )  =  Pu(r,  t )  (109) 

contains  only  high-wavenumber  information.  By  applying  V  to  Navier-Stokes  equa¬ 
tion  it  follows 

dU^’  —  +  u'(r,  t ) .  Vu'(r,  t )  -  V{p'/p]  +  y0V2u'(r,  t)  =  f(r,  t )  (110) 

where 

f(r,f)  =VF(r,t)  -T>[u(r,t)  •Vu(r,i)]  +  u'(r,t)-Vu'(r,t).  (Ill) 

Yuan  and  Ronis  (1992)  noted  that  this  ‘new’  force  contains  information  about 
boundaries  as  well  as  the  mode-coupling  effects  associated  with  velocity  components 
on  the  injection  scales.  Since  VF(r,  t)  =  0  away  from  boundaries,  they  defined  the 
random  stirring  force  used  in  RG  studies  to  result  from  the  mode  coupling  between  the 
energy  containing  and  the  smaller  scales.  Since  the  motion  on  all  scales  is  expected  to 
be  ‘chaotic’  including  the  energy  containing  scale,  Yuan  and  Ronis  (1992)  expected 
that  f(r,  t)  will  have  complicated,  chaotic  time  and  space  dependences.  Therefore, 
they  identified  this  as  the  quantity  which  is  actually  modeled  by  a  stochastic  force  in 
random  stirring  models  of  turbulence. 

The  autocorrelation  of  the  transverse  parts  of  f(r,  t)  is  assumed  with  a  non-zero 
correlation  time  (colored  noise).  The  major  problem  with  such  an  assumption  is  that 
the  resulting  theory  would  not  be  invariant  under  Galilean  transformation.  Yuan  and 
Ronis  (1992)  argued  that  there  is  no  a  priori  reason  why  global  Galilean  invariance 

must  hold.  Their  reason  is  that  the  random  forces  represents  the  effects  of  boundaries 
and  these  are  not  included  in  a  Galilean  transformation.  Note  that  Eq.  (110)  is  the 
subgrid  NSE  after  the  full  fluctuating  NSE  is  divided  into  the  super-  and  sub-  grid 
scales.  However,  the  filtering  operation  is  not  used  in  the  later  development,  and  in 
particular,  the  path-integral  RG  works  only  with  the  full  fluctuation  Navier-Stokes 
equation. 
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7.4  The  work  of  Carati 


The  objective  of  Carati  (1990a, b)  is  to  modify  the  e-expansion  of  YO. 

Carati  (1990a)  first  noticed  the  difference  between  the  e-  RG  applications  in  Ma 
and  Mazenko  (1975)  and  YO.  In  YO,  the  dimensionality  d  is  considered  as  a  fixed. 
A  stochastic  forcing  is  introduced  to  replace  the  initial  and  boundary  conditions  in 
Navier-Stokes  equation.  A  correlation  function  of  this  forcing  is  assumed  to  following 
a  powerlaw  with  exponent  y.  YO  then  defined  y  —  e  +  d  —  4  where  e-expansion  is 
performed  for  y.  Carati  (1990a)  attempted  to  modify  the  YO  procedure  by  introduc¬ 
ing  a  parametric  dimension  d  =  do  —  0(e).  But,  as  Carati  (1990a)  noted,  the  physical 
meaning  of  e?o  is  n°t  clear  a  prior. 

The  next  effort  of  Carati  (1990b)  involved  the  extension  of  the  white  noise  of 
the  forcing  correlation  function  to  a  colored  one  (again  a  colored  correlation  refers 
to  a  non-zero  correlation  in  time).  He  argued,  as  did  Yuan  and  Ronis  (1992),  that 
the  correlation  time  as  well  as  the  correlation  length  can  be  important  due  to  the 
appearance  of  macroscopic  scale  structures  in  turbulence.  It  is  then  possible  to  relate 
the  expansion  parameter  to  the  stochastic  forcing  correlation.  The  result  is  that  the 
parameter  e  is  not  treated  as  a  small  parameter  and  is  fixed  in  its  original  physical 
value  4.  The  frequency  correlation  provides  another  free  parameter  through  the  pow¬ 
erlaw  exponent  7/3  where  — 1  <  7  <  1.  The  original  YO’s  e-expansion  is  replaced  by 
arguing  that  the  frequency  integral  is  nearly  divergent  because  of  the  assumed  value 
of  7.  Alternatively,  the  e-expansion  has  now  been  replaced  by  a  7-expansion  where 
7  =  1  —  S.  This  d-expansion  was  interpreted  as  a  scheme  in  which  a  large  amount  of 
energy  is  injected  into  the  system  (Carati,  1990b). 

Although  Carati  viewed  the  colored  forcing  correlation  as  a  good  way  to  sidestep 
the  e-expansion,  he  did  not  consider  that  the  introduction  of  the  colored  noise  violates 
the  Galilean  invariance.  This  important  issue  was  pointed  out  by  Yuan  and  Ronis 
(1992)  in  their  analysis  (see  previous  subsection). 


7.5  The  work  of  Liang  and  Diamond 

Liang  and  Diamond  (1993)  conducted  a  careful  examination  on  the  applicability  of 

e-RG  to  two  dimensional  (2D)  flows.  They  examined  both  2D  fluid  and  2D  Magne- 
tohydrodynamicic  (MHD)  turbulence  but  our  discussion  will  be  limited  to  the  fluid 
turbulence  part  of  the  paper  only. 

The  most  important  difference  between  2D  and  3D  turbulence  is  the  direction  of 
energy  transfer.  There  is  a  dual-cascade  phenomena,  in  which  the  energy  is  trans- 
fered  from  large  to  small  scales,  while  enstrophy  (defined  as  mean-square  vorticity), 
is  transferred  from  large  scale  to  small  scales  (Kraichnan  and  Montgomery,  1980). 
Liang  and  Diamond  (1993)  found  that  when  a  proper  inertial  range  energy  spectrum 
is  assumed,  no  fixed  point  can  be  found  for  positive  or  negative  eddy  viscosity.  There- 
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fore,  e-RG  can  not  be  used  to  analyze  the  physics  of  2d  turbulence.  They  found  that 
the  reason  of  this  difficulty  can  be  traced  to  the  dual  energy  cascade. 

One  important  point  should  be  stressed.  Although  Liang  and  Diamond  (1993) 
followed  most  of  the  procedure  of  YO  e-RG,  they  did  not  apply  the  e-expansion. 
Instead,  they  used  the  proper  inertial  range  scaling  as  an  input. 

To  identify  the  reason  behind  the  failure  of  e-RG  in  2D,  Liang  and  Diamond  (1993) 
compared  it  to  the  EDQNM  closure.  In  e-RG  procedure,  all  calculations  are  one 
point  and  only  eddy  damping  term  appear.  The  nonlinear  noise  (nonlinear  coupling) 
is  replaced  by  adding  stochastic  forces.  The  major  problem  is  that  the  properties  of 
the  stochastic  forces  can  not  be  determined  recursively.  They  are  determined  instead 
by  the  properties  of  the  self-similarity  range  spectra,  which  in  reality,  are  the  goals 
of  the  calculation.  On  the  other  hand,  a  two-point  closure  calculation  determines  the 
form  of  the  spectrum.  Therefore,  both  the  eddy  damping  and  the  nonlinear  model 
coupling  terms  will  now  appear.  Furthermore,  while  the  e-RG  only  incorporate  the 
nonlocal  effects  of  the  small  scales  on  large  scales,  EDQNM  (and  DIA)  take  both 
effects  into  account,  in  principle. 

7.6  The  work  of  Lam 

Lam’s  approach  (Lam  1992),  although  belonging  to  the  FNS  school  of  activity,  is 
quite  different  from  the  other  members  of  that  group.  It  does  not  make  use  of  the 
A-expansion  procedure  used  by  FNS,  nor  the  correspondence  principle  and  the  e- 
expansion  procedure  used  in  the  YO  theory.  His  interpretation  of  e-RG  is  based  on 
phenomenological  approach.  Zhou  (1995)  has  shown  that  Lam’s  model  is  essentially 
the  physical  space  version  of  the  classical  closure  theory  (Leslie  and  Quarini,  1979) 
in  spectral  space  and  consider  the  corresponding  treatment  of  the  eddy  viscosity  and 
energy  backscatter. 


7.7  The  work  of  Eyink 

Eyink  (1994)  used  the  momentum-shell  RG  method  of  Kadanoff- Wilson  based  on  the 
Martin-Siggia-Rose  (1978)  field-theory  formulation  of  stochastic  dynamics. 

One  of  the  major  points  in  Eyink  (1994)  is  that,  contrary  to  the  claim  of  YO, 
the  higher-order  nonlinear  terms  generated  in  their  RG  analysis  are  not  irrelevant 
but  marginal  by  power-counting.  Because  of  this,  the  terms  neglected  by  YO  are 
not  necessarily  small  even  for  small  e.  The  problem  was  traced  to  issues  relevant  to 
Galilean-invariant  theories.  It  is  shown  not  to  occur  in  the  FNS  original  analysis  since 
they  were  dealing  with  fluctuation  dynamics  for  equilibrium  NS  fluids.  The  issue  of 
Galilean  invariance  has  also  been  discussed  by  Zhou  and  Vahala  (1993b). 

Eyink  (1994)  also  presented  his  objections  to  the  YO  applications  of  e-RG  to  tur¬ 
bulence  modeling.  First,  he  repeated  Kraichnan’s  (1987)  arguments  that  the  power- 
laws  derived  by  YO  can  be  simply  derived  from  dimensional  analysis.  Second,  he 
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stated,  as  Zhou  et  al.  (1988),  that  by  setting  e  =  4,  there  is  no  reason  to  believe 
that  additional  nonlinearities  are  negligible  or  insignificant  to  the  Physics.  Finally, 
the  YO  analysis  does  not  include  the  rescaling  procedure,  which  Eyink  (1994)  viewed 
as  vital  to  error  estimation.  In  fact,  Eyink  (1994)  even  stated  that  ‘YO  thory  is  not 
an  RG  analysis  at  all!”. 
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Part  IV 

RG-based  turbulence  modelling 

8  Modeling  of  Reynolds  stress  using  RG 


Turbulent  flows  of  practical  interest  involve  a  broad  spectrum  of  length  and  time 

scales  and  require  some  type  of  modeling  for  Reynolds  stresses.  The  most  widely 
used  are  the  two-equation  turbulence  models  based  on  the  transport  equations  for 
parameters  that  involve  the  length  and  the  time  scales  and  require  the  simplest  level 
of  Reynolds  stress  closure  that  does  not  depend  specifically  on  the  flow  geometry. 
In  its  standard  form  the  two-equation  Reynolds  stress  turbulence  models  involve  the 
turbulence  kinetic  energy  and  dissipation  based  on  a  Boussinesq  type  approximation 
involving  isotropic  eddy  viscosity.  Such  a  representation  of  turbulence  is  often  not 
effective  from  both  theoretical  as  well  as  phenomenological  point  of  view.  To  overcome 
this,  models  that  are  nonlinear  (i.e.,  quadratic)  in  the  mean  strain  rate  were  proposed 
in  the  form  of  a  constitutive  relation  (Speziale,  1987;  Yoshizawa,  1984).  These  models 
are  capable  of  predicting  the  anisotropy  in  the  Reynolds  stresses  in  complex  flows  but 
require  empirical  evaluation  of  the  model  constants  (Speziale,  1991). 

Two  studies  presented  here5  address  the  need  for  a  more  effective  approach  in  the 
development  of  two-equation  turbulence  models,  and  in  this  context  the  renormaliza¬ 
tion  group  (RG)  theory  based  models  are  considered  (for  a  summary  of  the  ICASE 
panel  discussion  on  RG  based  turbulence  modeling,  see  Zhou  and  Speziale  (1994)). 
These  models  fall  into  two  distinct  categories:  (a)  e  —  RG,  where  a  small  parameter 
e  is  introduced  into  the  exponent  of  the  forcing  correlation  function  (with  the  forcing 
function  being  introduced  into  the  momentum  equation),  and  the  theory  is  then  de¬ 
veloped  for  e  <  1  (Yakhot  and  Orszag,  1986);  and  (b)  recursion- RG,  which  does  not 
rely  on  an  e-expansion,  and  treats  explicitly  the  cubic  nonlinearities  introduced  into 
the  renormalized  momentum  equation  (Zhou,  et  al.  1988).  It  should  be  noted  that  in 
the  t  —  RG  while  all  constants  generated  are  evaluated  in  the  limit  e  <C  1,  at  the  same 
time,  all  exponents  that  are  e-dependent  are  evaluated  at  e  =  4.  In  fact,  e  =  4  is  re¬ 
quired  in  the  e—RG  to  recover  the  Kolmogorov  energy  spectrum  in  the  inertial  range. 
In  addition,  e  —  RG  theory  can  only  take  into  account  non-local  interactions  (Smith 
and  Reynolds,  1992).  On  the  other  hand,  recursion-AG  does  not  rely  on  an  e-expan- 
sion,  and  treats  explicitly  the  cubic  nonlinearities  introduced  into  the  renormalized 
momentum  equation.  Moreover,  recursion- AG  can  handle  both  local  and  non-local 
interactions.  Effects  such  as  the  cusp  behavior  in  the  transport  coefficients  (Zhou 

5Note  that  recently  Nagano  and  Itazu  (1997)  attempted  to  derived  an  eddy  viscosity  model  using 
the  iteractive  averaging  method 
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and  Vahala,  1993a)  which  are  due  to  local  interactions  and  the  cubic  nonlinearities 
are  recovered  in  these  theories.  As  in  e  -  AG,  the  eddy  viscosity  is  readily  deter¬ 
mined  from  the  solution  of  a  relatively  simple  differential  equation  (Zhou  and  Vahala, 
1993a).  However,  unlike  e  -  RG,  the  transport  coefficients  are  determined  over  the 
whole  resolvable  scales  and  not  just  in  the  wavenumber  limit  k  — >  0. 

The  contributions  to  the  Reynolds  stress  tensor  Ty  from  the  conventional  time¬ 
averaging  of  the  equations  of  motion  are  due  to  the  nonlinear  coupling  term  and  its 
interaction  with  the  velocity  field  (Zhou  et  al.,  1993).  In  particular,  the  Reynolds 
stress  tensor,  Ty  =  rj+  +  r~~,  where  the  ri++  -  part  arises  from  the  infrared  limit 
of  k  — *  0  and  is  due  to  the  uj  -  uj  distant  interaction  limit  while  the  Ty  -  part, 
arises  from  the  0  <  k  <  kc  spectrum  and  is  due  to  the  uj  —  Uj  local  interaction 
limit.  It  should  be  noted  that  in  e  -  RG  approach,  one  takes  the  large-scale  infrared 
limit  k  — >  0.  In  essence,  this  forces  a  spectral  gap  between  the  resolvable  part  of  the 
flow  field  and  the  small  unresolved  scales.  If  this  spectral  gap  were  somehow  present 
initially,  it  would  be  quickly  populated  in  just  a  few  eddy  turn  over  times.  Thus, 
retaining  only  the  distant  interactions  may  not  be  appropriate.  In  fact,  it  has  been 
shown  (Zhou  and  Vahala,  1993a)  that  the  energy  transfer  function  that  corresponds  to 
local  interactions  accounts  for  most  of  the  energy  flow  out  of  the  resolvable  scales.  It 
is,  therefore,  important  to  retain  both  local  and  nonlocal  interactions  in  the  modeling 
of  the  Reynolds  stress  and  this  can  be  readily  achieved  by  recursion-AG.  But  in  the 
e  -  RG  model,  the  Reynolds  stress  ry  =  rj+  and  is  obtained  purely  from  the  uj  -uj 
interaction  in  the  small  unresolved  scale  momentum  equation. 

Rubinstein  and  Barton  (1990)  have  derived  a  Reynolds  stress  model  using  the 
e  —  RG  method  (which  corresponds  to  the  infrared  limit  of  k,  — *  0)  of  the  following 
form: 

r++  =  -IkS13  +  vridjUi  +  dtUj)  -  ^T[CTl{daUidQUjy 

+Cr2{daUidjUa  +  dcUjdiU*)*  +  CraidiUcdjUc)*}.  (112) 

The  constants  CT i  =  0.034,  CT 2  =  0.104  and  CT 3  =  —0.014,  and  (. .  .)*  denotes 
the  deviatoric  part  of  the  expression  within  the  parenthesis.  The  first  two  terms 
correspond  to  the  linear  model  and  vT  =  C^K2 fe  is  the  isotropic  eddy- viscosity, 
where  e  is  the  turbulence  dissipation  and  w  0.09  based  on  empirical  data  from 
equilibrium  boundary  layer  flows.  The  above  model  is  quadratic  in  mean  strain  rate, 
includes  the  effect  of  convection  and  diffusion  and  is  qualitatively  similar  to  other 
second  order  models  (Speziale,  1991). 

Now,  if  one  follows  standard  recursion- AG  procedures  (Zhou  et  al.,  1993),  the 
relevant  part  of  the  small  scale  velocity  field  that  contributes  to  Ty  is  obtained 
in  the  wavenumber  space.  Transforming  back  to  the  physical  space  the  following 
algebraic  representation  is  obtained  (Zhou  et  al.,  1993): 

r+~(x)  =  Cm^-idaUidpU&Up  +i~j] 
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where  Ui  is  the  time-averaged  mean  velocity  (and  the  second  term  inside  the  square 
brackets  is  obtained  by  switching  the  indices).  However,  for  most  flows  of  interest 
(inhomogeneous  flows),  the  cut-off  wavenumber  kc  varies  considerably  throughout  the 
flow  domain.  Thus,  if  one  wishes  to  use  the  local  approximation,  as  well  as  retain 
kc  ~  0(1)  effects  so  as  to  obtain  an  algebraic  form  for  T”<(*),  then  the  coefficients 
CR 1  and  Cr2  will  be  functions  of  the  flow  quantities.  As  a  first  attempt  at  applying  this 
recursion- J?0  model  we  make  the  lowest  order  approximation  that  these  coefficients 
are  constants.  In  particular,  for  the  present  analysis,  Cr\  and  Cm  are  taken  to  be 
0.025  and  0.342  x  10-3,  respectively  (Zhou  et  al.,  1994). 

Combining  (112)  -  (113)  a  formal  expression  for  Reynolds  stress  which  includes 
both  the  local  and  nonlocal  interactions  may  be  obtained.  It  is  of  some  interest  to 
note  that  integrity  basis  representations  are  commonly  employed  to  represent  the 
anisotropic  part  of  the  Reynolds  stress  tensor  for  three  dimensional  turbulent  flows 
based  on  a  systematic  derivation  from  a  hierarchy  of  second-order  closure  models 
(Gatski  and  Speziale,  1993).  It  can  readily  be  shown  that  the  tensors  that  constitute 
the  integrity  basis  are  recovered  for  most  part  when  the  proposed  recursion- RG  model 
is  recast  appropriately  (Zhou  et  ah,  1994). 

The  above  expression  for  Reynolds  stress  (112)  -  (113)  are  to  be  used  along  with 
the  equations  of  motion  by  specifying  turbulent  kinetic  energy  and  dissipation.  In 
two-equation  turbulence  models,  this  closure  is  achieved  through  the  development  of 
transport  equations  for  the  turbulent  kinetic  energy  and  dissipation  -  quantities  that 
are  directly  related  to  the  length  and  time  scales  -  of  the  following  general  form: 

dtK  +  UjdjK  =  V-e  +  di[(v0  +  vT/aK)diK]  (114) 

dte  +  Ujdj£  =  CelV£/K  -  Ce2£2/K  +  diKvo  +  vT/ae)di£]  (115) 

where,  v  =  y0  +  vT  is  the  total  viscosity,  V  =  -Tij{dUi/dxj)  is  the  turbulence 
production,  £  is  the  scalar  turbulent  dissipation  rate.  The  quantities  CE i,  Ce2)  oik,  cx£ 
are  dimensionless  and  taken  to  be  1.44,  1.92,  1.0  and  1.3,  respectively,  consistent  with 
the  standard  form  of  the  two-equation  K-£  model  (based  on  empirical  data  obtained 
from  equilibrium  boundary  layer  flows).  We  note  that  with  various  approximations, 
the  £  —  RG-based  formulations  computed  these  constants  as  1.42,  1.68,  0.719  and 
0.719,  respectively  (Yakhot  and  Orszag,  1986). 

The  RG  theory  is  utilized  to  develop  Reynolds  stress  closure  models  for  the  pre¬ 
diction  of  turbulent  separated  flows.  The  combined  model  includes  both  the  local  and 
nonlocal  interaction  of  all  the  relevant  resolvable  scales.  The  ability  of  the  proposed 
model  to  accurately  predict  separated  flows  is  analyzed  from  a  combined  theoretical 
and  computational  standpoint  by  considering  turbulent  flow  past  a  backward  facing 
step  as  a  test  case.  The  results  obtained  based  on  detailed  computations  demonstrate 


that  the  RG  model  can  yield  very  good  predictions  for  the  turbulent  flow  of  an  in¬ 
compressible  viscous  fluid  over  a  backward-facing  step  (Zhou  et  al,  1994).  Thus,  in 
spite  of  its  well  known  deficiencies,  when  the  anisotropy  of  the  turbulent  stresses  are 
properly  accounted  for,  the  two-equation  turbulence  models  can  be  quite  effective  for 
the  prediction  of  turbulent  separated  flows. 

9  RG  based  K  —  e  model 

YO  derived  the  K  —  e  two  equation  model  using  e-RG  method.  Speziale  (1990)  found 
that  the  original  YO  model  performs  poorly  in  homogeneous  shear  flow.  The  value  Ct\ 
derived  in  YO  yields  excessively  large  growth  rate  for  the  turbulent  kinetic  energy  in 
homogeneous  shear  flow  in  comparision  to  both  physical  and  numerical  experiments 
(Speziale,  1991). 

Smith  and  Reynolds  (1992)  found  some  algebraic  error  in  the  original  derivation 
of  YO.  The  coefficients  of  the  dissipation  term  in  e  equation  is  not  in  good  agreement 
with  generally  accepted  values.  Furthermore,  YO’s  derivation  did  not  yield  a  term 
responsible  for  the  production  in  the  e  equation. 

The  original  derivation  of  YO  was  revised  by  Yakhot  and  Smith  (YS)  (1992)  by 
the  following  features: 

1.  The  ‘infrared  cutoff’  of  the  random  force,  <  f>f>  >=  0  when  0  <  k  <  Al 


<  fa(k,t)h(k',lS)  >=  D0k  yDa&{k)8{k+k')6{t-t'),  AL<k<cx>  =  0,  0  <  k  <  AL 

(116) 

This  property  is  needed  in  the  derivation  of  the  equation  for  the  mean  rate  of 
energy  dissipation  e  (YS). 

2.  The  input  of  energy  spectrum  for  the  interval  0  <  k  <  Al 

E(k)  ~  ka  (117) 

is  required  to  evaluate  the  integrals  (with  a  =  2). 

3.  Performing  a  Reynolds  decomposition  of  T\  =  —2vo (V/Ui)(V/uj)(Vz,Ui)  into 
mean  U  and  fluctuating  u  velocities. 

The  derivation  of  YO  and  YS  starts  from  dynamical  equations  for  the  homoge¬ 
neous  part  of  the  instantaneous  rate  of  energy  dissipation  per  unit  mass  e  =  vQ  (V/a*)2 

Ve  Ti 

^  \  /** 

—  +  UiVi£  =  2v0{W jUi) (Vj/i)  -  2y0 (V jU^VjU^iyiUi) 
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T2 

-  2  vl  (VyVjtti)8  — 2i/0  ( V j«i)  ( Vj  V  jp)  +  Vo ViVi£  (118) 

After  some  work,  at  the  stirred  fluids  at  the  long-time  and  large-distance  limit, 

the  e-RG  dissipation  equation  (YS,  1992)  is  found  to  be 

D&  =  Cei{e / KfajdjUi  -  Ce2e2/K  +  S^oa^e)  -  Tt  (119) 


where  C£ i  =  1.42,  C&  =  1.68  and 


71  =  2v0Sij 


.  dui  dui 


dxi  dxj 


(120) 


The  e  expansion  procedure  and  above  mentioned  assumptions  have  been  employed. 

The  above  e-RG  dissipation  equation  is  not  closed.  The  neglect  of  7?.  is  a  formally 
justified  approximation  at  high  Reynolds  number  if  the  hypothesis  of  local  isotropy 
is  invoked.  Durbin  and  Speziale  (1991)  have  questioned  the  validity  of  local  isotropy 
in  strongly  strained  turbulence  flows.  Yakhot  et  al.  (1992)  have  proposed  a  model 
where  7Z  =  7Z(rj),  where  the  standard  form  of  the  model  is  recovered  for  71  — >  0  in 
the  limit  of  weak  strains.  Note  that  Durbin  (1990)  has  already  developed  a  model 
for  the  production  of  dissipation  along  these  lines  that  was  quadratic  in  the  ratio 
of  production  to  the  dissipation  and,  hence,  quartic  in  77.  Lam  (1994)  published  a 
critique  on  the  YS  derivation. 

Iterating  the  expression  for  71  using  the  Navier-Stokes  equation  will  generate  a 
power  series 


OO  CTS 

7l  =  vTS3Y:rn(—r  (121) 

~rn  e 

n= 0 

where  S  =  (2 SijSij)l/2.  It  is  not  possible  to  evalute  the  summation  since  the  values 
of  coefficients  are  unknown. 

The  7Z  is  modeled  via  three  steps: 

1.  The  summation  is  performed  for  the  geometric  senes  for  every  three  terms. 

This  procedure  reduces  the  numbers  of  unknown  coefficients  to  one ,  (3. 

(122) 

2.  Assuming  that  the  fixed  point  value  r) o  =  4.38  of  the  homogeneous  shear  flow 

in  the  equilibrium  states,  is  invariant  to  dropping  all  terms  but  those  in  (122)  above, 
Yakhot  et  al.  (1992)  postulated  that 

(123) 

3.  One  now  further  assumes  that  the  isotropic  Reynolds  stress  tvj  —  —2CvKr)ij 
(Vij  =  SijK/s) 
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_  <W(  1  -  v/vo)  g3  _  g(l  -  g/go)  g  _  <7 
l+jftf  if  1+/V  if  iJ  y' 

The  undertermined  constant  0  =  0.012  for  the  von  Karman  constant  0.4. 
final  e-RG  dissipation  rate  transport  equation  is  given 


(124) 

The 


Dte  =  C^e/K^ijdjUi  -  C&£2/K  +  di{avTdie)  (125) 

where  the  coefficient  Cel  is  given  by 


^  -  v/vo) 

l TW~ 


(126) 


The  model  of  Yakhot  et  al.  (1992)  has  been  tested  for  homogeneous  shear  flows 
and  for  flow  over  a  backward  facing  step.  Excellent  results  are  obtained  in  both  cases. 
Recently,  Yakhot  and  Orszag  have  extended  the  model  and  applied  it  to  complex  flows 
using  the  FLUENT  code. 


Part  V 

Conclusion 


In  this  review,  we  explained  the  concepts  of  terms  renormalization  and  renormaliza¬ 
tion  group  by  referencing  to  various  physical  systems,  such  as  the  ising  model.  We 
then  present  a  comprehensive  review  on  applications  of  the  method  of  renormaliza¬ 
tion  group  to  turbulence.  These  parts  should  be  sufficient  for  readers  who  wishes  to 
get  a  balanced  view  of  RG  in  turbulence.  For  a  few  selected  approaches,  we  have 
provided  further  technical  details.  We  conclude  with  a  discussion  of  the  relevance 
and  application  of  renormalization  group  to  turbulence  modelling. 
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